streamline.f

Go to the documentation of this file.

!                   *****************
                    MODULE streamline
!                   *****************
!
!
!***********************************************************************
! BIEF 7.3
!***********************************************************************
!
!brief    MODULE FOR PARALLEL CHARACTERISTICS FOR TELEMAC2D (OR TELEMAC3D)
!
!note     JAJ PINXIT :: JACEK.JANKOWSKI@BAW.DE :: ENJOY!
!+          JAJ STARTED WELL FRI JUL  4 12:26:57 CEST 2008
!+          JAJ RUNNING WELL FRI JUL 11 16:20:21 CEST 2008
!+          JAJ RUNNING EVEN BETTER WED JUL 16 18:25:35 CEST 2008
!+          JAJ PRIVATISED MON JUL 21 10:41:04 CEST 2008
!+          JAJ OPTIMISED WED JUL 30 15:40:45 CEST 2008
!
!history  J-M HERVOUET (LNHE)
!+        31/07/2012
!+        V6P2
!+   Use of new array ELTCAR.
!+   Dimensions in SCARACT reviewed (previous confusion between
!+   NPOIN and NPLOT). INTENT(OUT) changed into INTENT(INOUT) in
!+   subroutine organise_char. NPOINT=NPLOT replaces NPOINT=NPOIN
!+   before the call to SCHAR41. ADD_CHAR11 and ADD_CHAR41 deleted.
!+   SCHAR11 and SCHAR41 simplified. Arguments removed in SCARACT.
!+   SCHAR12 and SCHAR13 added. DX,DY,DZ added to CHARAC_TYPE.
!+   More data saved when touching a solid boundary: XPLOT, YPLOT, DX, DY
!+   All this ensures strict equality of scalar and parallel runs !!!!
!+   However in 3D the vertical velocity is computed in Telemac-3D and
!+   has truncation errors, this will trigger differences.
!+   Now the sub-domain at the foot of the characteristic is returned
!+   if SCARACT called with argument POST=.TRUE.
!
!history  J-M HERVOUET (LNHE)
!+        08/01/2013
!+        V6P3
!+   Advection subroutines from Tomawac re-implemented here
!+   See SCHAR41_PER and SCHAR41_PER_4D
!+   A posteriori interpolation now possible in all cases.
!
!history  J-M HERVOUET (LNHE)
!+        28/01/2013
!+        V6P3
!+   Bug corrected in PREP_SENDBACK, IF(NCHARA.EQ.0) RETURN causes bugs
!+   when one processor has nothing to send, some arrays were not
!+   initialised.
!
!history  J-M HERVOUET (LNHE)
!+        22/02/2013
!+        V6P3
!+   Particle tracking in //. 3 subroutines added: send_particles,
!+   add_particle, del_particle, to be used by subroutine derive.
!
!history  J-M HERVOUET (LNHE)
!+        16/04/2013
!+        V6P3
!+   Case of successive uses of SCARACT. Bug corrected in the section
!+   calling organise_chars when NPLOT > LAST_NPLOT.
!
!history  J-M HERVOUET (LNHE)
!+        26/04/2013
!+        V6P3
!+   Organise_chars changed: new strategy of memory allocation: same
!+   maximum size for all processors, useful for some scenarios with
!+   particle. A sub-domain without initial particle must be able to
!+   receive one and needs memory for it.
!
!history  A. JOLY (EDF R&D, LNHE)
!+        22/05/2013
!+        V6P3
!+   Routines added to deal with the transport of algae. For 2D only.
!
!history  C. GOEURY (EDF R&D, LNHE)
!+        29/05/2013
!+        V7P0
!+   Routine SCHAR11_STO for stochastic diffusion in 2D.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        08/07/2014
!+        V7P0
!+   In subroutines SCHAR41_... updating of coordinates after touching a
!+   boundary was missing in RECVCHAR and caused a bug if the trajectory
!+   went to another subdomain just after. Look for 08/07/2014 to trace
!+   the corrections.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        14/08/2015
!+        V7P1
!+   Argument NRK added to SCARACT, it was hardcoded to 3. It is the
!+   number of sub-steps to do the characteristic pathline.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        31/08/2015
!+        V7P1
!+   Subroutine introduce_recvchar optimised (after results of a
!+   profiler saying we spent 12% of the time in it, but this seems
!+   wrong).
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        13/09/2016
!+        V7P2
!+   Checking more memory allocation, with more error messages and hints
!+   for solving problems.
!
!history  S.E.BOURBAN (HRW)
!+        21/03/2017
!+        V7P3
!+   Replacement of the DATA declarations by the PARAMETER associates
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_parallel
      USE interface_parallel, ONLY : p_max,p_sum,
     &  p_sync, p_mpi_type_free,p_mpi_alltoall,p_mpi_alltoallv
      IMPLICIT NONE
      PRIVATE
!
!     CALLED BY CHARAC (BIEF). CHARAC CALLED BY TELEMAC-2D AND 3D AND SISYPHE
!     CALLED BY THOMPS (TELEMAC-2D)
!     CALLED BY DERIVE (BIEF, USED BY TELEMAC-2D)
!     CALLED BY DERLAG (BIEF, USED BY TELEMAC-2D)
!     CALLED BY DERI3D (TELEMAC-3D)
!
      PUBLIC :: scaract,post_interp,send_particles,bief_interp
      PUBLIC :: add_particle,del_particle,organise_algs,send_info_alg
      PUBLIC :: del_info_alg
      PUBLIC :: oil_send_particles,oil_del_particle,oil_send_info
      PUBLIC :: oil_organise_chars
      PUBLIC :: dealloc_streamline
!
!     MAX_BASKET_SIZE IS THE NUMBER OF ADVECTED VARIABLES IN ONE PROCEDURE CALL
!
      INTEGER, PARAMETER :: max_basket_size=10 ! LARGE
!
!     SEE CALL GET_MPI_PARAMETERS IN SCARACT
!
      INTEGER last_nomb,last_nplot
!
!     THE TYPE FOR CHARACTERISTICS - LOST TRACEBACKS
!     DESCRIBES A TRACEBACK LEAVING A PARTITION TO ANOTHER ONE
!     FOR 2D WE USE 3D -> KNE AND ZP ARE OBSOLETE THEN
!
!     SEE ORG_CHARAC_TYPE1 IN LIBRARY PARALLEL
!     IT MUST BE THE SAME LOCAL TYPE
!
!     YA: NOW THE TYPE IS DECLARED IN DECLARATIONS_PARALLEL IN PARALLEL
!
!     TYPE CHARAC_TYPE
!       SEQUENCE   ! BUT SEEMS USELESS
!       INTEGER :: MYPID ! PARTITION OF THE TRACEBACK ORIGIN (HEAD)
!       INTEGER :: NEPID ! THE NEIGHBOUR PARTITION THE TRACEBACK ENTERS TO
!       INTEGER :: INE   ! THE LOCAL 2D ELEMENT NR THE TRACEBACK ENTERS IN THE NEIGBOUR PARTITION
!       INTEGER :: KNE   ! THE LOCAL LEVEL THE TRACEBACK ENTERS IN THE NEIGBOUR PARTITION
!       INTEGER :: IOR   ! THE POSITION OF THE TRAJECTORY -HEAD- IN MYPID [THE 2D/3D NODE OF ORIGIN]
!       INTEGER :: ISP   ! CURRENT RUNGE-KUTTA STEPS PASSED AS COLLECTED
!       INTEGER :: NSP   ! TOTAL RUNGE-KUTTA STEPS
!       INTEGER :: IFR   ! FREQUENCY
!       DOUBLE PRECISION :: XP,YP,ZP,FP  ! THE (X,Y,Z,F)-POSITION
!       DOUBLE PRECISION :: DX,DY,DZ,DF  ! THE CORRESPONDING DISPLACEMENTS
!       DOUBLE PRECISION :: BASKET(MAX_BASKET_SIZE) ! VARIABLES INTERPOLATED AT THE FOOT
!     END TYPE CHARAC_TYPE
!
!     THE CORRESPONDING MPI TYPE
!
      INTEGER :: characteristic
!
!     STRUCTURES FOR ALL-TO-ALL COMMUNICATION / SEND AND RECEIVE WITH COUNTERS
!     HEAP/SEND/RECVCOUNTS : COUNT THE NUMBER OF LOST TRACEBACKS PARTITION-WISE
!     S/RDISPLS : DISPLACEMENTS IN PARTITION-WISE SORTED SEND/RECVCHARS
!     HEAPCHAR : FOR SAVING INITIALLY LOST CHARACTERISTICS AND COLLECTING
!                THE IMPLANTED TRACEBACKS LOCALISED IN MY PARTITION
!     WHILE COLLECTING IS DONE IN HEAPCHARS, MOST ACTIVE OPERATIONS IN RECVCHAR
!     SENDCHAR REQUIRED DUE TO THE SPECIFIC SORTING FOR MPI_ALLTOALLV (OPTIMISE?)
!
      type(charac_type),ALLOCATABLE,DIMENSION(:)::heapchar
      type(charac_type),ALLOCATABLE,DIMENSION(:)::sendchar
      type(charac_type),ALLOCATABLE,DIMENSION(:)::recvchar
      INTEGER, ALLOCATABLE, DIMENSION(:)::sendcounts,sdispls
      INTEGER, ALLOCATABLE, DIMENSION(:)::recvcounts,rdispls
      INTEGER, ALLOCATABLE, DIMENSION(:)::heapcounts
!     STATIC DIMENSION FOR HEAPCHAR, SENDCHAR, RECVCHAR (SORRY, STATIC) in scaract
      INTEGER nchdim
!
!     WORK FIELD FOR COUNTING OCCURANCES PRO RANK / SORTING SENDCHAR
!
      INTEGER, ALLOCATABLE, DIMENSION(:) :: icha
!
!     STRUCTURE TO SEND THE INFO ASSOCIATED WITH ALGAE TRANSPORT
!
!
!     STRUCTURE TO SEND THE INFO ASSOCIATED WITH ALGAE TRANSPORT
!     THE STRUCTURE IS DECLARED IN DECLARATIONS_PARALLEL
!     YA: NOW THE TYPE IS DECLARED IN DECLARATIONS_PARALLEL IN PARALLEL
!
!     TYPE ALG_TYPE
!       SEQUENCE   ! NECESSARY TO DEFINE MPI TYPE ALG_CHAR
!       INTEGER :: MYPID ! PARTITION OF THE TRACEBACK ORIGIN (HEAD)
!       INTEGER :: NEPID ! THE NEIGHBOUR PARTITION THE TRACEBACK ENTERS TO
!       INTEGER :: IGLOB  ! THE GLOBAL NUMBER OF THE PARTICLES
!       INTEGER :: FLAG  ! USED TO ALIGN FIELDS
!       DOUBLE PRECISION :: VX,VY,VZ  ! THE (X,Y,Z) PARTICLE VELOCITY
!       DOUBLE PRECISION :: UX,UY,UZ  ! THE (X,Y,Z) FLUID VELOCITY
!       DOUBLE PRECISION :: UX_AV,UY_AV,UZ_AV  ! THE (X,Y,Z) AVERAGE FLUID VELOCITY
!       DOUBLE PRECISION :: K_AV,EPS_AV  ! THE VALUES OF K AND EPS
!       DOUBLE PRECISION :: H_FLU  ! THE WATER DEPTH AT POSITION OF VELOCITY
!       DOUBLE PRECISION :: PSI(3*101) ! VARIABLE PSI USED FOR THE BASSET FORCE
!     END TYPE ALG_TYPE
!
!     THE CORRESPONDING MPI TYPE
!
      INTEGER :: alg_char
!
!     STRUCTURES FOR ALL-TO-ALL COMMUNICATION / SEND AND RECEIVE WITH COUNTERS
!     HEAPALG : FOR SAVING INFORMATION TO BE SEND TO AT THE SAME TIME AS
!               PARTICLE POSITION. HEAP/SEND/RECVCOUNTS AND S/RDISPLS ARE USED
!               AS WELL. SENDALG AND RECVALG ARE USED IN A SIMILAR FASHION AS
!               SENDCHAR AND RECVCHAR
!
      TYPE(alg_type),ALLOCATABLE,DIMENSION(:)::heapalg
      TYPE(alg_type),ALLOCATABLE,DIMENSION(:)::sendalg
      TYPE(alg_type),ALLOCATABLE,DIMENSION(:)::recvalg
!
!     FOR OIL SPILLS
!
!     TYPE OIL_TYPE
!       SEQUENCE
!       INTEGER :: MYPID ! PARTITION OF THE TRACEBACK ORIGIN (HEAD)
!       INTEGER :: NEPID ! THE NEIGHBOUR PARTITION THE TRACEBACK ENTERS TO
!       INTEGER :: INE   ! THE LOCAL 2D ELEMENT NR THE TRACEBACK ENTERS IN THE NEIGBOUR PARTITION
!       INTEGER :: KNE   ! THE LOCAL LEVEL THE TRACEBACK ENTERS IN THE NEIGBOUR PARTITION
!       INTEGER :: IOR   ! THE POSITION OF THE TRAJECTORY -HEAD- IN MYPID [THE 2D/3D NODE OF ORIGIN]
!       INTEGER :: STATE   ! CURRENT RUNGE-KUTTA STEPS PASSED AS COLLECTED
!       INTEGER :: TPSECH   ! TOTAL RUNGE-KUTTA STEPS
!       INTEGER :: IFR   ! FREQUENCY
!       DOUBLE PRECISION :: SURFACE
!       DOUBLE PRECISION :: MASS0
!       DOUBLE PRECISION :: MASS
!       DOUBLE PRECISION :: MASS_EVAP
!       DOUBLE PRECISION :: MASS_DISS
!       DOUBLE PRECISION :: MASS_HAP(10)
!       DOUBLE PRECISION :: MASS_COMPO(10)
!       DOUBLE PRECISION :: TB_HAP(10)
!       DOUBLE PRECISION :: TB_COMPO(10)
!       DOUBLE PRECISION :: SOL_HAP(10)
!       DOUBLE PRECISION :: SOL_COMPO(10)
!     END TYPE OIL_TYPE
!
!     THE CORRESPONDING MPI TYPE
!
      INTEGER :: oil_charac
!
!     STRUCTURES FOR ALL-TO-ALL COMMUNICATION / SEND AND RECEIVE WITH COUNTERS
!     HEAP/SEND/RECVCOUNTS : COUNT THE NUMBER OF LOST TRACEBACKS PARTITION-WISE
!     S/RDISPLS : DISPLACEMENTS IN PARTITION-WISE SORTED SEND/RECVCHARS
!     HEAPCHAR : FOR SAVING INITIALLY LOST CHARACTERISTICS AND COLLECTING
!                THE IMPLANTED TRACEBACKS LOCALISED IN MY PARTITION
!     WHILE COLLECTING IS DONE IN HEAPCHARS, MOST ACTIVE OPERATIONS IN RECVCHAR
!     SENDCHAR REQUIRED DUE TO THE SPECIFIC SORTING FOR MPI_ALLTOALLV (OPTIMISE?)
!
      type(oil_type),ALLOCATABLE,DIMENSION(:)::heapoil
      type(oil_type),ALLOCATABLE,DIMENSION(:)::sendoil
      type(oil_type),ALLOCATABLE,DIMENSION(:)::recvoil
!
!     IF SET TO TRUE, EVERY DETAILED DEBUGGING IS SWITCHED ON
!
      LOGICAL :: trace=.false.
!
!     TO OPTIMISE PERIODICITY
!
      INTEGER :: i
      INTEGER, PARAMETER :: nplanmax=200
      INTEGER :: eta1(nplanmax) = (/ (i+1, i = 1, nplanmax) /)
!
!     Initialisation marker for SCARACT
!
      LOGICAL :: init=.true.
!
      CONTAINS
!
!---------------------------------------------------------------------
!   <<<<<<<<<<<<<<<<<< CHARACTERISTICS: PRIVATE >>>>>>>>>>>>>>>>>>
!---------------------------------------------------------------------
!
      SUBROUTINE deorg_charac_type
        INTEGER ier
        CALL p_mpi_type_free(characteristic,ier)
        RETURN
      END SUBROUTINE deorg_charac_type
!
      SUBROUTINE dealloc_streamline
        IF(ALLOCATED(heapchar)) DEALLOCATE(heapchar)
        IF(ALLOCATED(sendchar)) DEALLOCATE(sendchar)
        IF(ALLOCATED(recvchar)) DEALLOCATE(recvchar)
        IF(ALLOCATED(sendcounts)) DEALLOCATE(sendcounts)
        IF(ALLOCATED(sdispls)) DEALLOCATE(sdispls)
        IF(ALLOCATED(recvcounts)) DEALLOCATE(recvcounts)
        IF(ALLOCATED(rdispls)) DEALLOCATE(rdispls)
        IF(ALLOCATED(heapcounts)) DEALLOCATE(heapcounts)
        IF(ALLOCATED(icha)) DEALLOCATE(icha)
        IF(ALLOCATED(heapalg)) DEALLOCATE(heapalg)
        IF(ALLOCATED(sendalg)) DEALLOCATE(sendalg)
        IF(ALLOCATED(recvalg)) DEALLOCATE(recvalg)
        IF(ALLOCATED(heapoil)) DEALLOCATE(heapoil)
        IF(ALLOCATED(sendoil)) DEALLOCATE(sendoil)
        IF(ALLOCATED(recvoil)) DEALLOCATE(recvoil)
        init=.true.
        trace=.false.
        RETURN
      END SUBROUTINE dealloc_streamline
!
!---------------------------------------------------------------------
!     GET/SET ON DIMENSIONS AND COUNTERS
!---------------------------------------------------------------------
!
      INTEGER FUNCTION get_max_basket_size()
        get_max_basket_size = max_basket_size
      END FUNCTION get_max_basket_size
!
!---------------------------------------------------------------------
! ALLOCATE ALL STATIC FIELDS FOR ALL-TO-ALL COMMUNICATION
! PREPARE THE MPI_TYPE FOR LOST CHARACTERISTICS / TRACEBACKS
!---------------------------------------------------------------------
!
      SUBROUTINE organise_chars(NPARAM,NOMB,NCHDIM,LAST_NCHDIM) ! WATCH OUT
        USE bief_def, ONLY: ncsize
        IMPLICIT NONE
        INTEGER, INTENT(IN)    :: NPARAM,NOMB
        INTEGER, INTENT(INOUT) :: NCHDIM,LAST_NCHDIM
        INTEGER ISIZE
        IF (.NOT.ALLOCATED(heapcounts)) ALLOCATE(heapcounts(ncsize))
        IF (.NOT.ALLOCATED(sendcounts)) ALLOCATE(sendcounts(ncsize))
        IF (.NOT.ALLOCATED(recvcounts)) ALLOCATE(recvcounts(ncsize))
        IF (.NOT.ALLOCATED(sdispls))    ALLOCATE(sdispls(ncsize))
        IF (.NOT.ALLOCATED(rdispls))    ALLOCATE(rdispls(ncsize))
        IF (.NOT.ALLOCATED(icha))       ALLOCATE(icha(ncsize))
        heapcounts=0
        sendcounts=0
        recvcounts=0
        sdispls=0
        rdispls=0
        icha=0
!
        nchdim=nparam
!
        IF(last_nchdim.EQ.0) THEN
!         CASE OF A FIRST CALL (THINK THAT NCHDIM MAY BE 0)
          isize=max(nchdim,1)
          ALLOCATE(sendchar(isize))
          ALLOCATE(recvchar(isize))
          ALLOCATE(heapchar(isize))
!         LAST_NCHDIM WILL NEVER BE 0 AGAIN
          last_nchdim=isize
        ELSEIF(nchdim.GT.last_nchdim) THEN
!         NEW CALL WITH LARGER NCHDIM, WE HAVE TO REALLOCATE
          DEALLOCATE(sendchar)
          DEALLOCATE(recvchar)
          DEALLOCATE(heapchar)
          ALLOCATE(sendchar(nchdim))
          ALLOCATE(recvchar(nchdim))
          ALLOCATE(heapchar(nchdim))
          last_nchdim=nchdim
        ELSE
!         NEW CALL BUT DIMENSIONS OK, NOTHING TO DO
        ENDIF
        CALL org_charac_type1(nomb,characteristic) ! COMMIT THE CHARACTERISTICS TYPE FOR COMM
        RETURN
      END SUBROUTINE organise_chars
!
!---------------------------------------------------------------------
! ALLOCATE ALL STATIC FIELDS FOR TO FIND THE ELEMENT AND SUBDOMAIN
! IN ALGAE TRANSPORT
!---------------------------------------------------------------------
!
      SUBROUTINE organise_chars_for_a(NPARAM,NOMB,NCHDIM,LAST_NCHDIM)
        USE bief_def, ONLY: ncsize
        IMPLICIT NONE
        INTEGER, INTENT(IN)    :: NPARAM,NOMB
        INTEGER, INTENT(INOUT) :: NCHDIM,LAST_NCHDIM
        INTEGER ISIZE
!
!       IN CASE OF SEQUENTIAL CALCULATIONS
!
        isize=max(1,ncsize)
        IF (.NOT.ALLOCATED(heapcounts)) ALLOCATE(heapcounts(isize))
        IF (.NOT.ALLOCATED(sendcounts)) ALLOCATE(sendcounts(isize))
        IF (.NOT.ALLOCATED(recvcounts)) ALLOCATE(recvcounts(isize))
        IF (.NOT.ALLOCATED(sdispls))    ALLOCATE(sdispls(isize))
        IF (.NOT.ALLOCATED(rdispls))    ALLOCATE(rdispls(isize))
        IF (.NOT.ALLOCATED(icha))       ALLOCATE(icha(isize))
        heapcounts=0
        sendcounts=0
        recvcounts=0
        sdispls=0
        rdispls=0
        icha=0
!
        nchdim=nparam
!
        IF(last_nchdim.EQ.0) THEN
!         CASE OF A FIRST CALL (THINK THAT NCHDIM MAY BE 0)
          isize=max(nchdim,1)
          ALLOCATE(sendchar(isize))
          ALLOCATE(recvchar(isize))
          ALLOCATE(heapchar(isize))
!         LAST_NCHDIM WILL NEVER BE 0 AGAIN
          last_nchdim=isize
        ELSEIF(nchdim.GT.last_nchdim) THEN
!         NEW CALL WITH LARGER NCHDIM, WE HAVE TO REALLOCATE
          DEALLOCATE(sendchar)
          DEALLOCATE(recvchar)
          DEALLOCATE(heapchar)
          ALLOCATE(sendchar(nchdim))
          ALLOCATE(recvchar(nchdim))
          ALLOCATE(heapchar(nchdim))
          last_nchdim=nchdim
        ELSE
!         NEW CALL BUT DIMENSIONS OK, NOTHING TO DO
        ENDIF
!       TO ALLOW VARIABLES TO BE ALLOCATED IN SCALAR SIMULATIONS,
!       BUT WITHOUT USE OF MPI
        IF(ncsize.GT.1) THEN
          CALL org_charac_type1(nomb,characteristic) ! COMMIT THE CHARACTERISTICS TYPE FOR COMM
          ENDIF
        RETURN
      END SUBROUTINE organise_chars_for_a

!=================================================================================
!                            OILSPILL
!=================================================================================
!
      SUBROUTINE oil_organise_chars(NPARAM) ! WATCH OUT
          USE bief_def, ONLY: ncsize
          IMPLICIT NONE
          INTEGER, INTENT(IN)    :: NPARAM
          IF (.NOT.ALLOCATED(heapcounts)) ALLOCATE(heapcounts(ncsize))
          IF (.NOT.ALLOCATED(sendcounts)) ALLOCATE(sendcounts(ncsize))
          IF (.NOT.ALLOCATED(recvcounts)) ALLOCATE(recvcounts(ncsize))
          IF (.NOT.ALLOCATED(sdispls))    ALLOCATE(sdispls(ncsize))
          IF (.NOT.ALLOCATED(rdispls))    ALLOCATE(rdispls(ncsize))
          IF (.NOT.ALLOCATED(icha))       ALLOCATE(icha(ncsize))
          heapcounts=0
          sendcounts=0
          recvcounts=0
          sdispls=0
          rdispls=0
          icha=0
          ALLOCATE(sendoil(nparam))
          ALLOCATE(recvoil(nparam))
          ALLOCATE(heapoil(nparam))
          CALL org_charac_type_oil(oil_charac) ! COMMIT THE CHARACTERISTICS TYPE FOR COMM
          RETURN
        END SUBROUTINE oil_organise_chars

!=================================================================================
!                            OILSPILL
!=================================================================================

!
!---------------------------------------------------------------------
! ALLOCATE ALL STATIC FIELDS FOR ALL-TO-ALL COMMUNICATION
! PREPARE THE MPI_TYPE FOR ALGAE INFORMATION EXCHANGE
!---------------------------------------------------------------------
!
        SUBROUTINE organise_algs(NPARAM,NOMB)
          USE bief_def, ONLY: ncsize
          IMPLICIT NONE
          INTEGER, INTENT(IN)    :: NPARAM,NOMB
!
          IF (.NOT.ALLOCATED(heapcounts)) ALLOCATE(heapcounts(ncsize))
          IF (.NOT.ALLOCATED(sendcounts)) ALLOCATE(sendcounts(ncsize))
          IF (.NOT.ALLOCATED(recvcounts)) ALLOCATE(recvcounts(ncsize))
          IF (.NOT.ALLOCATED(sdispls))    ALLOCATE(sdispls(ncsize))
          IF (.NOT.ALLOCATED(rdispls))    ALLOCATE(rdispls(ncsize))
          IF (.NOT.ALLOCATED(icha))       ALLOCATE(icha(ncsize))
          heapcounts=0
          sendcounts=0
          recvcounts=0
          sdispls=0
          rdispls=0
          icha=0
          ALLOCATE(sendalg(nparam))
          ALLOCATE(recvalg(nparam))
          ALLOCATE(heapalg(nparam))
!         COMMIT THE CHARACTERISTICS TYPE FOR COMM
          CALL org_charac_type_alg(alg_char,nomb)
          RETURN
        END SUBROUTINE organise_algs
!
!---------------------------------------------------------------------
! FOR COLLECTING CHARACTERISTICS LEAVING INITIALLY A GIVEN PARTITION
! TO BE CALLED IN MODIFIED CHAR11 OR CHAR41 SIMILAR TO THE ORIGINAL
! BIEF SUBROUTINES / NOTE THE COUNTER NCHARA/HEAPCHAR USAGE
!---------------------------------------------------------------------
!
        SUBROUTINE collect_char(MYPID,IOR,MYII,IFACE,KNE,IFR,
     &                          ISP,NSP,XP,YP,ZP,FP,DX,DY,DZ,DF,
     &                          IFAPAR,NCHDIM,NCHARA)
          USE declarations_special
          IMPLICIT NONE
          INTEGER,  INTENT(IN) :: MYPID,IOR,MYII,IFACE,KNE,IFR
          INTEGER,  INTENT(IN) :: ISP,NSP,NCHDIM
          INTEGER,  INTENT(IN) :: IFAPAR(6,*)
          INTEGER,  INTENT(INOUT) :: NCHARA
          DOUBLE PRECISION, INTENT(IN) :: XP,YP,ZP,FP,DX,DY,DZ,DF
          INTEGER :: NEPID,II
          !
          IF(nchara.EQ.0) heapcounts=0
          nepid=ifapar(iface  ,myii)
          ii   =ifapar(iface+3,myii)
          nchara=nchara+1
          IF(nchara.GT.nchdim) THEN ! PROBABLY EXAGGERATED
            WRITE (lu,*) ' '
            WRITE (lu,*) 'MODULE STREAMLINE:'
            WRITE (lu,*) 'NCHARA=',nchara,' NCHDIM=',nchdim
            WRITE (lu,*) 'COLLECT_CHAR::NCHARA>NCHDIM, INCREASE'
            WRITE (lu,*) 'SECURITY COEFFICIENT FOR SCARACT'
            WRITE (lu,*) 'MYPID=',mypid
            CALL plante(1)
            stop
          ENDIF
          heapchar(nchara)%MYPID=mypid ! THE ORIGIN PID
          heapchar(nchara)%NEPID=nepid ! THE NEXT PID
          heapchar(nchara)%INE=ii      ! ELEMENT THERE
          heapchar(nchara)%KNE=kne     ! LEVEL THERE
          heapchar(nchara)%IOR=ior     ! THE ORIGIN 2D OR 3D NODE
          heapchar(nchara)%ISP=isp     ! R-K STEP AS COLLECTED
          heapchar(nchara)%NSP=nsp     ! R-K STEPS TO BE DONE
          heapchar(nchara)%IFR=ifr     ! FREQUENCY THERE
          heapchar(nchara)%XP=xp       ! X-POSITION
          heapchar(nchara)%YP=yp       ! Y-POSITION
          heapchar(nchara)%ZP=zp       ! Z-POSITION
          heapchar(nchara)%FP=fp       ! F-POSITION
          heapchar(nchara)%DX=dx       ! X-DISPLACEMENT
          heapchar(nchara)%DY=dy       ! Y-DISPLACEMENT
          heapchar(nchara)%DZ=dz       ! Z-DISPLACEMENT
          heapchar(nchara)%DF=df       ! F-DISPLACEMENT
!
          heapcounts(nepid+1)=heapcounts(nepid+1)+1
!
          RETURN
        END SUBROUTINE collect_char
!
!---------------------------------------------------------------------
! USED TO PUT THE ALGAE POSITION IN HEAPCHAR WHEN LOOKING FOR
! THE ELEMENT NUMBER AND PROCESSOR AFTER ALGAE TRANSPORT
!---------------------------------------------------------------------
!
        SUBROUTINE collect_alg(MYPID,NEPID,INE,KNE,ISP,NSP,
     &                         IFR,XP,YP,ZP,FP,DX,DY,DZ,DF,
     &                         NCHARA,NCHDIM)
          USE declarations_special
          IMPLICIT NONE
          INTEGER,  INTENT(IN) :: MYPID,NEPID,INE(*),KNE(*)
          INTEGER,  INTENT(IN) :: ISP,NSP,IFR,NCHARA,NCHDIM
          DOUBLE PRECISION, INTENT(IN) :: XP(*),YP(*),ZP(*),FP(*)
          DOUBLE PRECISION, INTENT(IN) :: DX(*),DY(*),DZ(*),DF(*)
          INTEGER I
!
          IF(nchara.GT.nchdim) THEN ! PROBABLY EXAGGERATED
            WRITE (lu,*) ' '
            WRITE (lu,*) 'MODULE STREAMLINE:'
            WRITE (lu,*) 'NCHARA=',nchara,' NCHDIM=',nchdim
            WRITE (lu,*) 'COLLECT_ALG::NCHARA>NCHDIM, INCREASE'
            WRITE (lu,*) 'SECURITY COEFFICIENT FOR SCARACT'
            WRITE (lu,*) 'MYPID=',mypid
            CALL plante(1)
            stop
          ENDIF
          IF(nchara.NE.0) THEN
            DO i=1,nchara
!             IN PARALLEL, ONLY HEAPCOUNTS WILL BE USED
!             EVEN IF RECVCOUNTS ARE DEFINED
!             NOTE: IN PARALLEL CALCULATIONS THE ORDER OF HEAPCHAR
!             NEEDS TO BE DEFINED BECAUSE OF THE WAY PREP_INITIAL_SEND
!             IS WRITTEN
              heapcounts(nepid+1)=i
              heapchar(i)%MYPID=mypid ! THE ORIGIN PID
              heapchar(i)%NEPID=nepid ! THE NEXT PID
              heapchar(i)%NEPID=nepid ! THE NEXT PID
              heapchar(i)%INE=ine(nchara-i+1)  ! ELEMENT THERE
              heapchar(i)%KNE=kne(nchara-i+1)  ! LEVEL THERE
              heapchar(i)%IOR=nchara-i+1       ! THE ORIGIN 2D OR 3D NODE
              heapchar(i)%ISP=isp     ! R-K STEP AS COLLECTED
              heapchar(i)%NSP=nsp     ! R-K STEPS TO BE DONE
              heapchar(i)%IFR=ifr     ! FREQUENCY THERE
              heapchar(i)%XP=xp(nchara-i+1)    ! X-POSITION
              heapchar(i)%YP=yp(nchara-i+1)    ! Y-POSITION
              heapchar(i)%ZP=zp(nchara-i+1)    ! Z-POSITION
              heapchar(i)%FP=fp(nchara-i+1)    ! F-POSITION
              heapchar(i)%DX=dx(nchara-i+1)    ! X-DISPLACEMENT
              heapchar(i)%DY=dy(nchara-i+1)    ! Y-DISPLACEMENT
              heapchar(i)%DZ=dz(nchara-i+1)    ! Z-DISPLACEMENT
              heapchar(i)%DF=df(nchara-i+1)    ! F-DISPLACEMENT
!             IN SCALAR MODE, ONLY RECVCOUNTS WILL BE USED
              recvcounts(nepid+1)=i
              recvchar(i)%MYPID=mypid ! THE ORIGIN PID
              recvchar(i)%NEPID=nepid ! THE NEXT PID
              recvchar(i)%INE=ine(i)  ! ELEMENT THERE
              recvchar(i)%KNE=kne(i)  ! LEVEL THERE
              recvchar(i)%IOR=i       ! THE ORIGIN 2D OR 3D NODE
              recvchar(i)%ISP=isp     ! R-K STEP AS COLLECTED
              recvchar(i)%NSP=nsp     ! R-K STEPS TO BE DONE
              recvchar(i)%IFR=ifr     ! FREQUENCY THERE
              recvchar(i)%XP=xp(i)    ! X-POSITION
              recvchar(i)%YP=yp(i)    ! Y-POSITION
              recvchar(i)%ZP=zp(i)    ! Z-POSITION
              recvchar(i)%FP=fp(i)    ! F-POSITION
              recvchar(i)%DX=dx(i)    ! X-DISPLACEMENT
              recvchar(i)%DY=dy(i)    ! Y-DISPLACEMENT
              recvchar(i)%DZ=dz(i)    ! Z-DISPLACEMENT
              recvchar(i)%DF=df(i)    ! F-DISPLACEMENT
            ENDDO
          ENDIF
!
          RETURN
        END SUBROUTINE collect_alg
!
!---------------------------------------------------------------------
! RE-INITIALISE THE STRUCTURE AFTER COMPLETING ALL ACTIONS
!---------------------------------------------------------------------
!
        SUBROUTINE re_initialise_chars(NSEND,NLOSTCHAR,NLOSTAGAIN,NARRV)
          IMPLICIT NONE
          INTEGER, INTENT(OUT) :: NSEND,NLOSTCHAR,NLOSTAGAIN,NARRV
          nlostchar=0
          nlostagain=0
          narrv=0
          nsend=0
!         MAYBE NOT MANDATORY
          IF (ALLOCATED(heapcounts)) heapcounts=0
          IF (ALLOCATED(sendcounts)) sendcounts=0
          IF (ALLOCATED(recvcounts)) recvcounts=0
          IF (ALLOCATED(sdispls))    sdispls=0    ! NOT NECESSARY?
          IF (ALLOCATED(rdispls))    rdispls=0    ! NOT NECESSARY?
          IF (ALLOCATED(icha))      icha=0        ! NOT NECESSARY?
        END SUBROUTINE re_initialise_chars
!
!---------------------------------------------------------------------
! PREPARE THE INITIAL SEND OF THE LOST CHARACTERISTICS
! THE FIELDS ARE PREPARED ACCORDING THE MPI_ALLTOALL(V) REQUIREMENTS
!---------------------------------------------------------------------
!
        SUBROUTINE prep_initial_send(NSEND,NLOSTCHAR,NCHARA)
          USE bief_def, ONLY : ncsize
          IMPLICIT NONE
          INTEGER, INTENT(IN)    :: NSEND
          INTEGER, INTENT(OUT)   :: NLOSTCHAR
          INTEGER, INTENT(INOUT) :: NCHARA
          INTEGER I,N
!         TODO: THIS LINE IS PERHAPS A BUG, SEE PREP_SENDBACK...
          IF(nchara.EQ.0) THEN
            nlostchar = 0
            RETURN
          ENDIF
          sendcounts=heapcounts
          sdispls(1) = 0 ! CONTIGUOUS DATA
          DO i=2,ncsize
            sdispls(i) = sdispls(i-1)+sendcounts(i-1)
          ENDDO
          icha=sendcounts ! A RUNNING COUNTER PARTITION-WISE
          DO i=1,nchara
!           HEAPCHAR(I)%NEPID+1 - THE PARTITION WE SEND TO / OR -1
            IF(heapchar(i)%NEPID.GE.0) THEN
              n=heapchar(i)%NEPID+1
              sendchar(sdispls(n)+icha(n))=heapchar(i)
              icha(n)=icha(n)-1
            ENDIF
          ENDDO
          nlostchar = nsend
          heapcounts=0
          nchara=0
          RETURN
        END SUBROUTINE prep_initial_send
!
!---------------------------------------------------------------------
! PREPARE HEAPALG AND SENDALG ACCORDING TO THE MPI_ALLTOALL(V)
! REQUIREMENTS
!---------------------------------------------------------------------
!
        SUBROUTINE prep_initial_send_alg(NSEND,NLOSTCHAR,NCHARA)
          USE bief_def, ONLY : ncsize
          IMPLICIT NONE
          INTEGER, INTENT(IN)    :: NSEND
          INTEGER, INTENT(OUT)   :: NLOSTCHAR
          INTEGER, INTENT(INOUT) :: NCHARA
          INTEGER I,N
!         TODO: THIS LINE IS PERHAPS A BUG, SEE PREP_SENDBACK...
          IF(nchara.EQ.0) RETURN
          sendcounts=heapcounts
          sdispls(1) = 0 ! CONTIGUOUS DATA
          DO i=2,ncsize
            sdispls(i) = sdispls(i-1)+sendcounts(i-1)
          END DO
          icha=sendcounts ! A RUNNING COUNTER PARTITION-WISE
          DO i=1,nchara
!           HEAPCHAR(I)%NEPID+1 - THE PARTITION WE SEND TO / OR -1
            IF(heapalg(i)%NEPID.GE.0) THEN
              n=heapalg(i)%NEPID+1
              sendalg(sdispls(n)+icha(n))=heapalg(i)
              icha(n)=icha(n)-1
            ENDIF
          ENDDO
          nlostchar = nsend
          heapcounts=0
          nchara=0
          RETURN
        END SUBROUTINE prep_initial_send_alg

!====================================================================
!                         OILSPILL
!====================================================================

        SUBROUTINE oil_prep_initial_send(NSEND,NLOSTCHAR,NCHARA)
          USE bief_def, ONLY : ncsize
          IMPLICIT NONE
          INTEGER, INTENT(IN)    :: NSEND
          INTEGER, INTENT(OUT)   :: NLOSTCHAR
          INTEGER, INTENT(INOUT) :: NCHARA
          INTEGER I,N
!         TODO: THIS LINE IS PERHAPS A BUG, SEE PREP_SENDBACK...
          IF(nchara.EQ.0) RETURN
          sendcounts=heapcounts
          sdispls(1) = 0 ! CONTIGUOUS DATA
          DO i=2,ncsize
            sdispls(i) = sdispls(i-1)+sendcounts(i-1)
          END DO
          icha=sendcounts ! A RUNNING COUNTER PARTITION-WISE
          DO i=1,nchara
!           HEAPCHAR(I)%NEPID+1 - THE PARTITION WE SEND TO / OR -1
            IF(heapoil(i)%NEPID.GE.0) THEN
              n=heapoil(i)%NEPID+1
              sendoil(sdispls(n)+icha(n))=heapoil(i)
              icha(n)=icha(n)-1
            ENDIF
          ENDDO
          nlostchar = nsend
          heapcounts=0
          nchara=0
          RETURN
        END SUBROUTINE oil_prep_initial_send

!====================================================================
!                          OILSPILL
!====================================================================

!
!---------------------------------------------------------------------
! COLLECT IMPLANTED TRACEBACKS WHICH ARE COMPLETED/LOCALISED
! ON A HEAP, SETTING BY THE WAY ALSO THE NUMBER OF THE LOST-AGAIN
! TRACEBACKS ACCORDINGLY TO THE PARTITION THEY SHOULD BE SEND TO
! THE "LOCALISED" MARK IS SET IN ADD_CHAR11/41
!---------------------------------------------------------------------
!
        SUBROUTINE heap_found(NLOSTAGAIN,NARRV,NCHARA)
          USE declarations_special
          IMPLICIT NONE
          INTEGER, INTENT(OUT)    :: NLOSTAGAIN
          INTEGER, INTENT(IN)     :: NARRV
          INTEGER, INTENT(INOUT)  :: NCHARA
          INTEGER I
          sendcounts=0
!         DO NOT ZEROIZE NCHARA, HEAPCOUNTS / ADDING FROM GENERATIONS!
!         COUNTER PARTITION-WISE, ALSO MY-OWN
          DO i=1,narrv
            IF(recvchar(i)%NEPID.EQ.-1) THEN ! A LOCALISED TRACEBACK
              nchara=nchara+1
              IF(nchara.GT.nchdim) THEN
                WRITE (lu,*) ' '
                WRITE(lu,*) 'MODULE STREAMLINE:'
                WRITE(lu,*) 'NOT ENOUGH MEMORY ALLOCATION, INCREASE'
                WRITE(lu,*) 'SECURITY COEFFICIENT FOR SCARACT'
                CALL plante(1)
                stop
              ENDIF
              heapchar(nchara) = recvchar(i) ! ALREADY INTERPOLATED?
              heapcounts(heapchar(nchara)%MYPID+1) =
     &             heapcounts(heapchar(nchara)%MYPID+1)+1
            ELSE ! A LOST-AGAIN CHARACTERISTIC / TO BE SORTED LATER
              sendcounts(recvchar(i)%NEPID+1) =
     &           sendcounts(recvchar(i)%NEPID+1)+1
            ENDIF
          END DO
          nlostagain=sum(sendcounts)
          RETURN
        END SUBROUTINE heap_found
!
!---------------------------------------------------------------------
! PREPARE LOST-AGAIN TRACEBACKS FOR THE NEXT COMMUNICATION
! FILL IN THE STRUCTURE FOR THE ALL-TO-ALL COMMUNICATION
! NOTE THAT SENDCOUNTS ARE SET IN HEAP_FOUND
! (OPTIMISE(?): HEAP-FOUND AND PREP_LOST_AGAIN CAN BE JOINED)
!---------------------------------------------------------------------
!
        SUBROUTINE prep_lost_again(NSEND,NARRV)
          USE bief_def, ONLY : ncsize
          IMPLICIT NONE
          INTEGER,INTENT(IN)   :: NARRV
          INTEGER, INTENT(OUT) :: NSEND
          INTEGER I,N
          sdispls(1) = 0 ! CONTIGUOUS DATA MARKER
          DO i=2,ncsize
            sdispls(i) = sdispls(i-1)+sendcounts(i-1)
          ENDDO
          icha=0
          nsend=0
          DO i=1,narrv
            n=recvchar(i)%NEPID
            IF (n.NE.-1) THEN  ! A LOST-AGAIN TRACEBACK
              icha(n+1)=icha(n+1)+1
              nsend=nsend+1
              sendchar(sdispls(n+1)+icha(n+1)) = recvchar(i)
            ENDIF
          ENDDO
          RETURN
        END SUBROUTINE prep_lost_again
!
!---------------------------------------------------------------------
! MOVE THE HEAP OF IMPLANTED AND COMPLETED TRACEBACKS (I.E. COMPLETED
! IN MY PARTITION) TO DATA STRUCTURES FOR ALL-TO-ALL COMMUNICATION
! SENDCHAR IS FILLED ACCORDING TO THE MPI_ALLTOALLV REQUIREMENTS
! ALL-TO-ALL PATTERN INCLUDE MY OWN LOST TRACEBACKS THAT CAME BACK
!---------------------------------------------------------------------
!
        SUBROUTINE prep_sendback(NCHARA)
          USE bief_def, ONLY: ipid, ncsize
          IMPLICIT NONE
          INTEGER, INTENT(INOUT) :: NCHARA
          INTEGER  :: I,N
!
          sendcounts=heapcounts
          sdispls(1) = 0 ! CONTIGUOUS DATA
          DO i=2,ncsize
            sdispls(i) = sdispls(i-1)+sendcounts(i-1)
          ENDDO
          icha=0
          IF(nchara.GT.0) THEN
            DO i=1,nchara
!             MYPID+1 - IS THE -ORIGIN- PARTITION
              n=heapchar(i)%MYPID+1
              icha(n)=icha(n)+1
              sendchar(sdispls(n)+icha(n))=heapchar(i)
!             SIGN IN THE SENDBACK ORIGIN FOR DEBUGGING PURPOSES
              sendchar(sdispls(n)+icha(n))%NEPID=ipid
            ENDDO
          ENDIF
          heapcounts=0
          nchara=0
          RETURN
        END SUBROUTINE prep_sendback
!
!---------------------------------------------------------------------
! THE GLOBAL COMMUNICATION OF LOST CHARACTERISTICS - ALL-TO-ALL
! (THIS IS THE HEART OF ALL THINGS / THE GLOBAL COMMUNICATION)
! THE DATA IS SENT AND (NOTE!) RECEIVED -SORTED- ACCORDING TO THE
! MPI_ALLTOALL(V) SPECIFICATION IN A CONTIGUOUS FIELDS
! DATA FOR A GIVEN PROCESSOR/PARTITION IN FIELD SECTIONS DESCRIBED BY
! DISPLACEMENTS SDISPLS AND RDISPLS
!---------------------------------------------------------------------
!
        SUBROUTINE glob_char_comm()
          USE bief_def, ONLY : ncsize
          USE declarations_special
          IMPLICIT NONE
          INTEGER :: I,IER
          CALL p_mpi_alltoall(sendcounts,1,mpi_integer,
     &                        recvcounts,1,mpi_integer,
     &                        ier)
          IF(ier.NE.mpi_success) THEN
            WRITE(lu,*)
     &       'STREAMLINE::GLOB_CHAR_COMM::MPI_ALLTOALL ERROR: ',ier
            CALL plante(1)
            stop
          ENDIF
          rdispls(1) = 0 ! SAVE THE RECEIVED DATA CONTIGUOUSLY
          DO i=2,ncsize
            rdispls(i) = rdispls(i-1)+recvcounts(i-1)
          ENDDO
          CALL p_mpi_alltoallv
     &      (sendchar,sendcounts,sdispls,characteristic,
     &       recvchar,recvcounts,rdispls,characteristic,
     &       ier)
          IF(ier.NE.mpi_success) THEN
            WRITE(lu,*)
     &       'STREAMLINE::GLOB_CHAR_COMM::MPI_ALLTOALLV ERROR: ',ier
            CALL plante(1)
            stop
          ENDIF
          RETURN
        END SUBROUTINE glob_char_comm
!
!---------------------------------------------------------------------
! THE GLOBAL COMMUNICATION OF ALGAE INFO - ALL-TO-ALL
! THE DATA IS SENT AND (NOTE!) RECEIVED -SORTED- ACCORDING TO THE
! MPI_ALLTOALL(V) SPECIFICATION IN A CONTIGUOUS FIELDS
! DATA FOR A GIVEN PROCESSOR/PARTITION IN FIELD SECTIONS DESCRIBED BY
! DISPLACEMENTS SDISPLS AND RDISPLS
!---------------------------------------------------------------------
!
        SUBROUTINE glob_alg_comm()
          USE bief_def, ONLY : ncsize
          USE declarations_special
          IMPLICIT NONE
          INTEGER :: I,IER
!
          CALL p_mpi_alltoall(sendcounts,1,mpi_integer,
     &                        recvcounts,1,mpi_integer,
     &                        ier)
!
          IF(ier.NE.mpi_success) THEN
            WRITE(lu,*)
     &       ' @STREAMLINE::GLOB_CHAR_COMM::MPI_ALLTOALL ERROR: ',ier
            CALL plante(1)
            stop
          ENDIF
          rdispls(1) = 0 ! SAVE THE RECEIVED DATA CONTIGUOUSLY
          DO i=2,ncsize
            rdispls(i) = rdispls(i-1)+recvcounts(i-1)
          ENDDO
          CALL p_mpi_alltoallv
     &      (sendalg,sendcounts,sdispls,alg_char,
     &       recvalg,recvcounts,rdispls,alg_char,
     &       ier)

          IF(ier.NE.mpi_success) THEN
            WRITE(lu,*)
     &       ' @STREAMLINE::GLOB_ALG_COMM::MPI_ALLTOALLV ERROR: ',ier
            CALL plante(1)
            stop
          ENDIF
!
          RETURN
        END SUBROUTINE glob_alg_comm

!====================================================================
!                         OILSPILL
!====================================================================

        SUBROUTINE oil_glob_char_comm()
          USE bief_def, ONLY : ncsize
          USE declarations_special
          IMPLICIT NONE
          INTEGER :: I,IER
          CALL p_mpi_alltoall(sendcounts,1,mpi_integer,
     &                        recvcounts,1,mpi_integer,
     &                        ier)
          IF(ier.NE.mpi_success) THEN
            WRITE(lu,*)
     &       ' @STREAMLINE::GLOB_CHAR_COMM::MPI_ALLTOALL ERROR: ',ier
            CALL plante(1)
            stop
          ENDIF
          rdispls(1) = 0 ! SAVE THE RECEIVED DATA CONTIGUOUSLY
          DO i=2,ncsize
            rdispls(i) = rdispls(i-1)+recvcounts(i-1)
          ENDDO
          CALL p_mpi_alltoallv
     &      (sendoil,sendcounts,sdispls,oil_charac,
     &       recvoil,recvcounts,rdispls,oil_charac,
     &       ier)
          IF(ier.NE.mpi_success) THEN
            WRITE(lu,*)
     &       ' @STREAMLINE::GLOB_CHAR_COMM::MPI_ALLTOALLV ERROR: ',ier
            CALL plante(1)
            stop
          ENDIF
          RETURN
        END SUBROUTINE oil_glob_char_comm
!
!====================================================================
!                         OILSPILL
!====================================================================
!

!---------------------------------------------------------------------
! TELEMAC3D PRISMS, INTERPOLATION OF RECVCHAR
! ELT,ETA AND SHP,SHZ MUST BE CORRECTLY PROVIDED VIA ADD_CHAR11
!   -> MATCHED TO THE RANGE 1:NRANGE (NO CHECKING - FOR SPEED!!!)
!   N IS THE POSITION FORESEEN FOR A GIVEN VARIABLE VAL IN THE BASKET
!---------------------------------------------------------------------
!
        SUBROUTINE interp_recvchar_41
     &    (val,n,ikle,elt,eta,fre,shp,shz,shf,nelem,npoin2,
     &     nplan,nrange,post,nomb,perio,ya4d)
          USE declarations_special
          IMPLICIT NONE
          INTEGER, INTENT(IN) :: N,NELEM,NPOIN2,NPLAN,NRANGE,NOMB
          INTEGER, INTENT(IN) :: IKLE(nelem,3)
          INTEGER, INTENT(IN) :: ELT(nrange),ETA(nrange)
          INTEGER, INTENT(IN) :: FRE(nrange)
          DOUBLE PRECISION, INTENT(IN) :: SHP(3,nrange),SHZ(nrange)
          DOUBLE PRECISION, INTENT(IN) :: SHF(*)
          DOUBLE PRECISION, INTENT(IN) :: VAL(npoin2,nplan,*)
          LOGICAL, INTENT(IN) :: POST,PERIO,YA4D
          INTEGER I,ETAP1,I1,I2,I3,IFR
          DOUBLE PRECISION UMSHZ,UMSHF
!
!         INTERPOLATION
!
          IF(nomb.GT.0) THEN
            IF(perio) THEN
              eta1(nplan)=1
              IF(ya4d) THEN
                DO i=1,nrange
                  IF(recvchar(i)%NEPID.EQ.-1) THEN ! LOCALISED
                    i1=ikle(elt(i),1)
                    i2=ikle(elt(i),2)
                    i3=ikle(elt(i),3)
                    umshz=1.d0-shz(i)
                    umshf=1.d0-shf(i)
                    etap1=eta1(eta(i))
                    ifr=fre(i)
                    recvchar(i)%BASKET(n) = umshf *
     &        ((val(i1,eta(i),ifr  ) * shp(1,i)
     &        + val(i2,eta(i),ifr  ) * shp(2,i)
     &        + val(i3,eta(i),ifr  ) * shp(3,i)) * umshz
     &       +( val(i1,etap1 ,ifr  ) * shp(1,i)
     &        + val(i2,etap1 ,ifr  ) * shp(2,i)
     &        + val(i3,etap1 ,ifr  ) * shp(3,i)) * shz(i) )
     &                 + shf(i) *
     &        ((val(i1,eta(i),ifr+1) * shp(1,i)
     &        + val(i2,eta(i),ifr+1) * shp(2,i)
     &        + val(i3,eta(i),ifr+1) * shp(3,i)) * umshz
     &       +( val(i1,etap1 ,ifr+1) * shp(1,i)
     &        + val(i2,etap1 ,ifr+1) * shp(2,i)
     &        + val(i3,etap1 ,ifr+1) * shp(3,i)) * shz(i) )
!                 THIS IS JUST TO AVOID A BUG ON HP COMPILER
                  ELSEIF(recvchar(i)%NEPID.LT.-1) THEN
                    WRITE(lu,*) 'STREAMLINE INTERP_RECVCHAR_11'
                    WRITE(lu,*) 'NEPID OUT OF RANGE:',recvchar(i)%NEPID
                    WRITE(lu,*) 'FOR I=',i
                    CALL plante(1)
                    stop
!                 END OF THIS IS JUST TO AVOID A BUG ON HP COMPILER
                  ENDIF
                ENDDO
              ELSE
                DO i=1,nrange
                  IF(recvchar(i)%NEPID.EQ.-1) THEN ! LOCALISED
                    etap1=eta1(eta(i))
                    recvchar(i)%BASKET(n) =
     &          (val(ikle(elt(i),1),eta(i),1)*shp(1,i)
     &          +val(ikle(elt(i),2),eta(i),1)*shp(2,i)
     &          +val(ikle(elt(i),3),eta(i),1)*shp(3,i))*(1.d0-shz(i))
     &         +(val(ikle(elt(i),1),etap1 ,1)*shp(1,i)
     &          +val(ikle(elt(i),2),etap1 ,1)*shp(2,i)
     &          +val(ikle(elt(i),3),etap1 ,1)*shp(3,i))*shz(i)
!                 THIS IS JUST TO AVOID A BUG ON HP COMPILER
                  ELSEIF(recvchar(i)%NEPID.LT.-1) THEN
                    WRITE(lu,*) 'STREAMLINE INTERP_RECVCHAR_11'
                    WRITE(lu,*) 'NEPID OUT OF RANGE:',recvchar(i)%NEPID
                    WRITE(lu,*) 'FOR I=',i
                    CALL plante(1)
                    stop
!                 END OF THIS IS JUST TO AVOID A BUG ON HP COMPILER
                  ENDIF
                ENDDO
              ENDIF
!             RESTORING THE ORIGINAL ETA1
              eta1(nplan)=nplan+1
            ELSE
              DO i=1,nrange
              IF(recvchar(i)%NEPID.EQ.-1) THEN ! LOCALISED
                recvchar(i)%BASKET(n) =
     &          ( val(ikle(elt(i),1),eta(i),1)  *shp(1,i)
     &          + val(ikle(elt(i),2),eta(i),1)  *shp(2,i)
     &          + val(ikle(elt(i),3),eta(i),1)  *shp(3,i))*(1.d0-shz(i))
     &         +( val(ikle(elt(i),1),eta(i)+1,1)*shp(1,i)
     &          + val(ikle(elt(i),2),eta(i)+1,1)*shp(2,i)
     &          + val(ikle(elt(i),3),eta(i)+1,1)*shp(3,i))*shz(i)
!             THIS IS JUST TO AVOID A BUG ON HP COMPILER
              ELSEIF(recvchar(i)%NEPID.LT.-1) THEN
                WRITE(lu,*) 'STREAMLINE  INTERP_RECVCHAR_11'
                WRITE(lu,*) 'NEPID OUT OF RANGE:',recvchar(i)%NEPID
                WRITE(lu,*) 'FOR I=',i
                CALL plante(1)
                stop
!             END OF THIS IS JUST TO AVOID A BUG ON HP COMPILER
              ENDIF
              ENDDO
            ENDIF
          ENDIF
!
!         SAVING INTERPOLATION DATA
!
          IF(post) THEN
            IF(ya4d) THEN
              DO i=1,nrange
                IF(recvchar(i)%NEPID.EQ.-1) THEN ! LOCALISED
                  recvchar(i)%XP=shp(1,i)
                  recvchar(i)%YP=shp(2,i)
                  recvchar(i)%ZP=shp(3,i)
                  recvchar(i)%DX=shz(i)
                  recvchar(i)%DY=shf(i)
                  recvchar(i)%INE=elt(i)
                  recvchar(i)%KNE=eta(i)
                  recvchar(i)%IFR=fre(i)
                ENDIF
              ENDDO
            ELSE
              DO i=1,nrange
                IF(recvchar(i)%NEPID.EQ.-1) THEN ! LOCALISED
                  recvchar(i)%XP=shp(1,i)
                  recvchar(i)%YP=shp(2,i)
                  recvchar(i)%ZP=shp(3,i)
                  recvchar(i)%DX=shz(i)
                  recvchar(i)%INE=elt(i)
                  recvchar(i)%KNE=eta(i)
                ENDIF
              ENDDO
            ENDIF
          ENDIF
!
          RETURN
        END SUBROUTINE interp_recvchar_41
!
!---------------------------------------------------------------------
! TELEMAC2D TRIANGLES, INTERPOLATION OF RECVCHAR
! ELT AND SHP MUST BE CORRECTLY PROVIDED VIA SCHAR11
!   -> MATCHED TO THE RANGE 1:NRANGE (NO CHECKING - FOR SPEED!!!)
!   N IS THE POSITION FORESEEN FOR A GIVEN VARIABLE VAL IN THE BASKET
!---------------------------------------------------------------------
!
        SUBROUTINE interp_recvchar_11
     &      (val,n,ikle,elt,shp,nelem,npoin,nrange,ielm,post,nomb)
          USE bief
          USE declarations_special
          IMPLICIT NONE
          INTEGER, INTENT(IN) :: N,NELEM,NPOIN,NRANGE,IELM,NOMB
          INTEGER, INTENT(IN) :: IKLE(nelem,*)
          INTEGER, INTENT(IN) :: ELT(nrange)
          DOUBLE PRECISION, INTENT(IN) :: SHP(3,nrange)
          DOUBLE PRECISION, INTENT(IN) :: VAL(npoin)
          LOGICAL, INTENT(IN) :: POST
!
          DOUBLE PRECISION SHP11,SHP12,SHP14
          DOUBLE PRECISION SHP22,SHP23,SHP24
          DOUBLE PRECISION SHP33,SHP31,SHP34
!
!         SHOULD BE SAME EPSILO THAN SCHAR11 AND INTERP
          DOUBLE PRECISION, PARAMETER :: EPSILO = 1.d-6
          INTEGER I
!
          IF(nomb.GT.0) THEN
          IF(ielm.EQ.11) THEN
            DO i=1,nrange
              IF(recvchar(i)%NEPID.EQ.-1) THEN ! LOCALISED
                recvchar(i)%BASKET(n) =
     &                         val(ikle(elt(i),1)) * shp(1,i)
     &                       + val(ikle(elt(i),2)) * shp(2,i)
     &                       + val(ikle(elt(i),3)) * shp(3,i)
!             THIS IS JUST TO AVOID A BUG ON HP COMPILER
              ELSEIF(recvchar(i)%NEPID.LT.-1) THEN
                WRITE(lu,*) 'STREAMLINE INTERP_RECVCHAR_11'
                WRITE(lu,*) 'NEPID OUT OF RANGE:',recvchar(i)%NEPID
                WRITE(lu,*) 'FOR I=',i
                CALL plante(1)
                stop
!             END OF THIS IS JUST TO AVOID A BUG ON HP COMPILER
              ENDIF
            ENDDO
          ELSEIF(ielm.EQ.12) THEN
            DO i=1,nrange
              IF(recvchar(i)%NEPID.EQ.-1) THEN ! LOCALISED
                shp11=shp(1,i)-shp(3,i)
                shp12=shp(2,i)-shp(3,i)
                shp14=3.d0*shp(3,i)
                shp22=shp(2,i)-shp(1,i)
                shp23=shp(3,i)-shp(1,i)
                shp24=3.d0*shp(1,i)
                shp33=shp(3,i)-shp(2,i)
                shp31=shp(1,i)-shp(2,i)
                shp34=3.d0*shp(2,i)
!               SEE INTERP
                IF(     shp11.GT.    -2.d0*epsilo .AND.
     &                  shp11.LT.1.d0+4.d0*epsilo .AND.
     &                  shp12.GT.    -2.d0*epsilo .AND.
     &                  shp12.LT.1.d0+4.d0*epsilo .AND.
     &                  shp14.LT.1.d0+4.d0*epsilo ) THEN
                        recvchar(i)%BASKET(n) =
     &                      val(ikle(elt(i),1)) * shp11
     &                    + val(ikle(elt(i),2)) * shp12
     &                    + val(ikle(elt(i),4)) * shp14
                ELSEIF( shp22.GT.    -2.d0*epsilo .AND.
     &                  shp22.LT.1.d0+4.d0*epsilo .AND.
     &                  shp23.GT.    -2.d0*epsilo .AND.
     &                  shp23.LT.1.d0+4.d0*epsilo .AND.
     &                  shp24.LT.1.d0+4.d0*epsilo ) THEN
                        recvchar(i)%BASKET(n) =
     &                      val(ikle(elt(i),2)) * shp22
     &                    + val(ikle(elt(i),3)) * shp23
     &                    + val(ikle(elt(i),4)) * shp24
                ELSEIF( shp33.GT.    -2.d0*epsilo .AND.
     &                  shp33.LT.1.d0+4.d0*epsilo .AND.
     &                  shp31.GT.    -2.d0*epsilo .AND.
     &                  shp31.LT.1.d0+4.d0*epsilo .AND.
     &                  shp34.LT.1.d0+4.d0*epsilo ) THEN
                        recvchar(i)%BASKET(n) =
     &                      val(ikle(elt(i),3)) * shp33
     &                    + val(ikle(elt(i),1)) * shp31
     &                    + val(ikle(elt(i),4)) * shp34
                ELSE
                  WRITE(lu,*) 'I=',i,' NOT LOCATED, ELT=',elt(i)
                  CALL plante(1)
                  stop
                ENDIF
!             THIS IS JUST TO AVOID A BUG ON HP COMPILER
              ELSEIF(recvchar(i)%NEPID.LT.-1) THEN
                WRITE(lu,*) 'STREAMLINE  INTERP_RECVCHAR_11'
                WRITE(lu,*) 'NEPID OUT OF RANGE:',recvchar(i)%NEPID
                WRITE(lu,*) 'FOR I=',i
                CALL plante(1)
                stop
!             END OF THIS IS JUST TO AVOID A BUG ON HP COMPILER
              ENDIF
            ENDDO
          ELSEIF(ielm.EQ.13) THEN
          DO i=1,nrange
            IF(recvchar(i)%NEPID.EQ.-1) THEN ! LOCALISED
              recvchar(i)%BASKET(n) =
     &      val(ikle(elt(i),1)) * (2.d0*shp(1,i)-1.d0)* shp(1,i)
     &     +val(ikle(elt(i),2)) * (2.d0*shp(2,i)-1.d0)* shp(2,i)
     &     +val(ikle(elt(i),3)) * (2.d0*shp(3,i)-1.d0)* shp(3,i)
     &     +val(ikle(elt(i),4)) * 4.d0 * shp(1,i)*shp(2,i)
     &     +val(ikle(elt(i),5)) * 4.d0 * shp(2,i)*shp(3,i)
     &     +val(ikle(elt(i),6)) * 4.d0 * shp(3,i)*shp(1,i)
!           THIS IS JUST TO AVOID A BUG ON HP COMPILER
            ELSEIF(recvchar(i)%NEPID.LT.-1) THEN
              WRITE(lu,*) 'STREAMLINE  INTERP_RECVCHAR_11'
              WRITE(lu,*) 'NEPID OUT OF RANGE:',recvchar(i)%NEPID
              WRITE(lu,*) 'FOR I=',i
              CALL plante(1)
              stop
!           END OF THIS IS JUST TO AVOID A BUG ON HP COMPILER
            ENDIF
          ENDDO
          ELSE
            WRITE(lu,*) 'WRONG IELM IN INTERP_RECVCHAR_11:',ielm
            CALL plante(1)
            stop
          ENDIF
          ENDIF
!
          IF(post) THEN
            IF(ielm.EQ.11.OR.ielm.EQ.12.OR.ielm.EQ.13) THEN
              DO i=1,nrange
                IF(recvchar(i)%NEPID.EQ.-1) THEN ! LOCALISED
                  recvchar(i)%XP=shp(1,i)
                  recvchar(i)%YP=shp(2,i)
                  recvchar(i)%ZP=shp(3,i)
                  recvchar(i)%INE=elt(i)
                ENDIF
              ENDDO
            ELSE
              WRITE(lu,*) 'WRONG IELM IN INTERP_RECVCHAR_11:',ielm
              CALL plante(1)
              stop
            ENDIF
          ENDIF
!
          RETURN
        END SUBROUTINE interp_recvchar_11
!
!---------------------------------------------------------------------
! INTRODUCE RECEIVED VALUES IN THE BASKET BACK IN THE TELEMAC
! STRUCTURES, N BEING THE POSITION IN THE BASKET FOR A GIVEN VARIABLE
!---------------------------------------------------------------------
!
        SUBROUTINE introduce_recvchar(VAL,NOMB,NARRV,IELM,
     &                                SHP,SHZ,SHF,ELT,ETA,FRE,
     &                                POST,YA4D)
          USE bief
          IMPLICIT NONE
          INTEGER, INTENT(IN)    :: NOMB,NARRV,IELM
          INTEGER, INTENT(INOUT) :: ELT(*),ETA(*),FRE(*)
          DOUBLE PRECISION, INTENT(INOUT) :: SHP(3,*),SHZ(*)
          DOUBLE PRECISION, INTENT(INOUT) :: SHF(*)
!         IF YES, SAVE INTERPOLATION DATA
          LOGICAL, INTENT(IN) :: POST,YA4D
          TYPE(bief_obj), INTENT(INOUT) :: VAL
          INTEGER I,N,IPOIN,MAXDIM
          IF(nomb.GT.max_basket_size) THEN
            WRITE(lu,*) 'STREAMLINE::INTRODUCE_RECVCHAR'
            WRITE(lu,*) 'NOMB>MAX_BASKET_SIZE'
            WRITE(lu,*) 'NOMB,MAX_BASKET_SIZE=',nomb,max_basket_size
            CALL plante(1)
            stop
          ENDIF
!
!         SAVING INTERPOLATION DATA, DEPENDING ON IELM
!
          IF(post) THEN
            IF(ielm.EQ.11) THEN
              DO i=1,narrv
                ipoin=recvchar(i)%IOR
                shp(1,ipoin)=recvchar(i)%XP
                shp(2,ipoin)=recvchar(i)%YP
                shp(3,ipoin)=recvchar(i)%ZP
                elt(ipoin)  =recvchar(i)%INE
              ENDDO
            ELSEIF(ielm.EQ.41) THEN
              IF(ya4d) THEN
                DO i=1,narrv
                  ipoin=recvchar(i)%IOR
                  shp(1,ipoin)=recvchar(i)%XP
                  shp(2,ipoin)=recvchar(i)%YP
                  shp(3,ipoin)=recvchar(i)%ZP
                  shz(ipoin)  =recvchar(i)%DX
                  shf(ipoin)  =recvchar(i)%DY
                  elt(ipoin)  =recvchar(i)%INE
                  eta(ipoin)  =recvchar(i)%KNE
                  fre(ipoin)  =recvchar(i)%IFR
                ENDDO
              ELSE
                DO i=1,narrv
                  ipoin=recvchar(i)%IOR
                  shp(1,ipoin)=recvchar(i)%XP
                  shp(2,ipoin)=recvchar(i)%YP
                  shp(3,ipoin)=recvchar(i)%ZP
                  shz(ipoin)  =recvchar(i)%DX
                  elt(ipoin)  =recvchar(i)%INE
                  eta(ipoin)  =recvchar(i)%KNE
                ENDDO
              ENDIF
            ELSE
              WRITE(lu,*)'STREAMLINE::INTRODUCE_RECVCHAR'
              WRITE(lu,*)'UNEXPECTED IELM=',ielm
              CALL plante(1)
              stop
            ENDIF
          ENDIF
!
!         NOW THE INTERPOLATION, DEPENDING ON TYPE
!
          IF(nomb.GT.0) THEN
            IF(val%TYPE.EQ.2) THEN
              DO i=1,narrv
                val%R(recvchar(i)%IOR)=recvchar(i)%BASKET(1)
              ENDDO
            ELSEIF(val%TYPE.EQ.4) THEN
!             LOOKING FOR THE LARGEST SIZE IN ALL ARRAYS TO INTERPOLATE
              maxdim=val%ADR(1)%P%DIM1
              IF(nomb.GT.1) THEN
                DO n=2,nomb
                  maxdim=max(maxdim,val%ADR(n)%P%DIM1)
                ENDDO
              ENDIF
              DO n=1,nomb
                IF(maxdim.GT.val%ADR(n)%P%DIM1) THEN
!                 HERE LOSS OF OPTIMISATION:
!                 IF CHARACTERISTICS COMPUTED FOR QUASI BUBBLE
!                 OR QUADRATIC, FUNCTIONS WHICH ARE LINEAR
!                 MUST NOT BE INTERPOLATED BEYOND THEIR SIZE.
                  DO i=1,narrv
                    IF(recvchar(i)%IOR.LE.val%ADR(n)%P%DIM1) THEN
                      val%ADR(n)%P%R(recvchar(i)%IOR)=
     &                recvchar(i)%BASKET(n)
                    ENDIF
                  ENDDO
                ELSE
!                 HERE NO RISK
                  DO i=1,narrv
                    val%ADR(n)%P%R(recvchar(i)%IOR)=
     &              recvchar(i)%BASKET(n)
                  ENDDO
                ENDIF
              ENDDO
            ELSE
              WRITE(lu,*)'STREAMLINE::INTRODUCE_RECVCHAR'
              WRITE(lu,*)'UNEXPECTED VAL%TYPE=',val%TYPE
              CALL plante(1)
              stop
            ENDIF
          ENDIF
          RETURN
        END SUBROUTINE introduce_recvchar

!/////////////////////////////////////////////////////////////////////
!
!            TELEMAC ROUTINES MODIFIED FOR THE PURPOSE OF
!           PARALLEL STREAMLINE TRACKING - MOSTLY FROM BIEF
!
!/////////////////////////////////////////////////////////////////////

!                       ******************
                        SUBROUTINE schar41
!                       ******************
!
     &( u , v , w , dt , nrk , x , y , zstar , z , ikle2 , ibor ,
     &  xplot , yplot , zplot , dx , dy , dz , shp , shz , elt , eta ,
     &  nplot , npoin2 , nelem2 ,nelmax2, nplan , surdet ,
     &  sens  , ifapar, nchdim,nchara,add,sigma)
!
!***********************************************************************
! BIEF VERSION 7.3
!
! 08/11/04 : ADAPTATION A LA TRANSFORMEE SIGMA GENERALISEE
! 12/10/05 : BUG CORRIGE, VOIR VARIABLE IELE QUI ETAIT AVANT IEL
!            ET EFFACAIT UN AUTRE IEL.
!
!
!***********************************************************************
!
!  FONCTION :
!
!     REMONTEE OU DESCENTE
!     DES COURBES CARACTERISTIQUES
!     SUR DES PRISMES DE TELEMAC-3D
!     DANS L'INTERVALLE DE TEMPS DT
!     AVEC UNE DISCRETISATION ELEMENTS FINIS
!
!
!  DISCRETISATION :
!
!     LE DOMAINE EST APPROCHE PAR UNE DISCRETISATION ELEMENTS FINIS
!     UNE APPROXIMATION LOCALE EST DEFINIE POUR LE VECTEUR VITESSE :
!     LA VALEUR EN UN POINT D'UN ELEMENT NE DEPEND QUE DES VALEURS
!     AUX NOEUDS DE CET ELEMENT
!
!
!  RESTRICTIONS ET HYPOTHESES :
!
!     LE CHAMP CONVECTEUR U EST SUPPOSE INDEPENDANT DU TEMPS
!
!-----------------------------------------------------------------------
!                             ARGUMENTS
! .________________.____.______________________________________________.
! |      NOM       |MODE|                   ROLE                       |
! |________________|____|______________________________________________|
! |    U,V,W       | -->| COMPONENTS OF ADVECTION VELOCITY
! |                |    | BUT W IS W* x DELTAZ (STEMS FROM TRIDW2)
! |    DT          | -->| PAS DE TEMPS.                                |
! |    NRK         | -->| NOMBRE DE SOUS-PAS DE RUNGE-KUTTA.           |
! |    X,Y,ZSTAR   | -->| COORDONNEES DES POINTS DU MAILLAGE.          |
! |    Z           | -->| COTE DANS LE MAILLAGE REEL                   |
! |    IKLE2       | -->| TRANSITION ENTRE LES NUMEROTATIONS LOCALE    |
! |                |    | ET GLOBALE DU MAILLAGE 2D.                   |
! |    IBOR        | -->| NUMEROS 2D DES ELEMENTS AYANT UNE FACE COMMUNE
! |                |    | AVEC L'ELEMENT .  SI IFABOR<0 OU NUL         |
! |                |    | ON A UNE FACE LIQUIDE,SOLIDE,OU PERIODIQUE   |
! |  X..,Y..,ZPLOT |<-->| POSITIONS SUCCESSIVES DES DERIVANTS.         |
! |    DX,DY,DZ    | -- | STOCKAGE DES SOUS-PAS .                      |
! |    SHP         |<-->| COORDONNEES BARYCENTRIQUES 2D AU PIED DES    |
! |                |    | COURBES CARACTERISTIQUES.                    |
! |    SHZ         |<-->| COORDONNEES BARYCENTRIQUES SUIVANT Z DES     |
! |                |    | NOEUDS DANS LEURS ETAGES "ETA" ASSOCIES.     |
! |    ELT         |<-->| NUMEROS DES ELEMENTS 2D CHOISIS POUR CHAQUE  |
! |                |    | NOEUD.                                       |
! |    ETA         |<-->| NUMEROS DES ETAGES CHOISIS POUR CHAQUE NOEUD.|
! |    NPLOT       | -->| NOMBRE DE DERIVANTS.                         |
! |    NPOIN2      | -->| NOMBRE DE POINTS DU MAILLAGE 2D.             |
! |    NELEM2      | -->| NOMBRE D'ELEMENTS DU MAILLAGE 2D.            |
! |    NPLAN       | -->| NOMBRE DE PLANS.
! |    SIGMA       | -->| IF YES, TRANSFORMED MESH
! |    SURDET      | -->| VARIABLE UTILISEE PAR LA TRANSFORMEE ISOPARAM.
! |    SENS        | -->| DESCENTE OU REMONTEE DES CARACTERISTIQUES.   |
! |    ISO         | -- | STOCKAGE BINAIRE DE LA FACE DE SORTIE.       |
! |________________|____|______________________________________________|
!  MODE: -->(DONNEE NON MODIFIEE),<--(RESULTAT),<-->(DONNEE MODIFIEE)
!-----------------------------------------------------------------------
!     - APPELE PAR : CARACT
!     - PROGRAMMES APPELES : NEANT
!
!***********************************************************************
!
      USE bief
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)    :: SENS,NPLAN,NCHDIM,NELMAX2
      INTEGER         , INTENT(IN)    :: NPOIN2,NELEM2,NPLOT,NRK
      INTEGER         , INTENT(IN)    :: IKLE2(nelmax2,3)
      INTEGER         , INTENT(INOUT) :: ELT(nplot),NCHARA
      DOUBLE PRECISION, INTENT(IN)    :: U(npoin2,nplan),V(npoin2,nplan)
      DOUBLE PRECISION, INTENT(IN)    :: W(npoin2,nplan),SURDET(nelem2)
      DOUBLE PRECISION, INTENT(INOUT) :: XPLOT(nplot),YPLOT(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: ZPLOT(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: SHP(3,nplot),SHZ(nplot)
      DOUBLE PRECISION, INTENT(IN)    :: X(npoin2),Y(npoin2),DT
      DOUBLE PRECISION, INTENT(IN)    :: Z(npoin2,nplan),ZSTAR(nplan)
      DOUBLE PRECISION, INTENT(INOUT) :: DX(nplot),DY(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: DZ(nplot)
      INTEGER         , INTENT(IN)    :: IBOR(nelmax2,5,nplan-1)
      INTEGER         , INTENT(INOUT) :: ETA(nplot)
      INTEGER         , INTENT(IN)    :: IFAPAR(6,*)
      LOGICAL, INTENT(IN)             :: ADD,SIGMA
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELE,ISO,ISPDONE,NSP
      INTEGER :: IPLOT,ISP,I1,I2,I3,IEL,IET,IET2,ISOH,ISOV,IFA,ISUI(3)
!
      DOUBLE PRECISION PAS,A1,DX1,DY1,DXP,DYP,XP,YP,ZP,NUM,DENOM
      DOUBLE PRECISION DELTAZ,PAS2,ZUP,ZDOWN,ZZ
!
      INTRINSIC abs , int , max , sqrt
!
      parameter( isui = (/ 2 , 3 , 1 /) )
      DOUBLE PRECISION, PARAMETER :: EPSILO = -1.d-6
      DOUBLE PRECISION, PARAMETER :: EPSDZ = 1.d-4
!
!-----------------------------------------------------------------------
!     FOR EVERY POINT
!-----------------------------------------------------------------------
!
      DO iplot = 1 , nplot
!
        IF(add) THEN
!
        xplot(iplot)   = recvchar(iplot)%XP
        yplot(iplot)   = recvchar(iplot)%YP
        zplot(iplot)   = recvchar(iplot)%ZP
        dx(iplot)      = recvchar(iplot)%DX
        dy(iplot)      = recvchar(iplot)%DY
        dz(iplot)      = recvchar(iplot)%DZ
        elt(iplot)     = recvchar(iplot)%INE
        eta(iplot)     = recvchar(iplot)%KNE
        nsp            = recvchar(iplot)%NSP ! R-K STEPS TO BE FULLFILLED
        ispdone        = recvchar(iplot)%ISP ! R-K STEPS ALREADY DONE
        iel = elt(iplot)
        iet = eta(iplot)
        xp  = xplot(iplot)
        yp  = yplot(iplot)
        zp  = zplot(iplot)
        i1 = ikle2(iel,1)
        i2 = ikle2(iel,2)
        i3 = ikle2(iel,3)
        shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                 -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
        shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                 -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
        shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                 -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
        IF(sigma) THEN
          shz(iplot) = (zp-zstar(iet)) / (zstar(iet+1)-zstar(iet))
        ELSE
          zdown=shp(1,iplot)*z(i1,iet)
     &         +shp(2,iplot)*z(i2,iet)
     &         +shp(3,iplot)*z(i3,iet)
          zup  =shp(1,iplot)*z(i1,iet+1)
     &         +shp(2,iplot)*z(i2,iet+1)
     &         +shp(3,iplot)*z(i3,iet+1)
          shz(iplot) = (zp-zdown) / max(zup-zdown,epsdz)
        ENDIF
!       ASSUME ALL ARE LOCALISED, IT WILL BE SET OTHERWISE IF LOST-AGAIN
        recvchar(iplot)%NEPID=-1
!
        ELSE
!
          iel = elt(iplot)
!         POINTS WITH IEL=0 ARE TREATED SO THAT THE FINAL
!         INTERPOLATION GIVES 0.,
!         AND WE SKIP TO NEXT POINT IPLOT (CYCLE)
          IF(iel.EQ.0) THEN
            elt(iplot)=1
            eta(iplot)=1
            shp(1,iplot)=0.d0
            shp(2,iplot)=0.d0
            shp(3,iplot)=0.d0
            shz(iplot)=0.d0
            cycle
          ENDIF
          iet = eta(iplot)
          i1 = ikle2(iel,1)
          i2 = ikle2(iel,2)
          i3 = ikle2(iel,3)
          dxp =( u(i1,iet  )*shp(1,iplot)
     &         + u(i2,iet  )*shp(2,iplot)
     &         + u(i3,iet  )*shp(3,iplot) )*(1.d0-shz(iplot))
     &        +( u(i1,iet+1)*shp(1,iplot)
     &         + u(i2,iet+1)*shp(2,iplot)
     &         + u(i3,iet+1)*shp(3,iplot) )*shz(iplot)
          dyp =( v(i1,iet  )*shp(1,iplot)
     &         + v(i2,iet  )*shp(2,iplot)
     &         + v(i3,iet  )*shp(3,iplot) )*(1.d0-shz(iplot))
     &        +( v(i1,iet+1)*shp(1,iplot)
     &         + v(i2,iet+1)*shp(2,iplot)
     &         + v(i3,iet+1)*shp(3,iplot) )*shz(iplot)
!         VERTICAL VELOCITY NOT CONSIDERED HERE !!
          nsp=max(1,int(nrk*dt*sqrt((dxp**2+dyp**2)*surdet(iel))))
          ispdone=1
!
        ENDIF
!
        pas = sens * dt / nsp
!
!       LOOP ON RUNGE-KUTTA SUB-STEPS
!
!       COMPILER MUST DO NOTHING IF ISPDONE>NSP
!       IN MODE "ADD", ISP = ISPDONE HAS NOT BEEN FULLY DONE
!       IT IS RESTARTED HERE
!
      DO isp = ispdone,nsp
!
!-----------------------------------------------------------------------
!       LOCALISING THE ARRIVAL POINT
!-----------------------------------------------------------------------
!
        pas2 = pas
!
!       IN MODE "ADD" ITERATIONS ALREADY DONE ARE SKIPPED AND
!                     CHARACTERISTICS GONE IN ANOTHER SUB-DOMAIN SKIPPED
!
        IF(add) THEN
          IF(isp.EQ.ispdone) GO TO 50
          IF(recvchar(iplot)%NEPID.NE.-1) cycle
        ENDIF
!
        iel = elt(iplot)
        iet = eta(iplot)
        i1 = ikle2(iel,1)
        i2 = ikle2(iel,2)
        i3 = ikle2(iel,3)
!
        dx(iplot) = ((u(i1,iet  )*shp(1,iplot)
     &              + u(i2,iet  )*shp(2,iplot)
     &              + u(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &              +(u(i1,iet+1)*shp(1,iplot)
     &              + u(i2,iet+1)*shp(2,iplot)
     &              + u(i3,iet+1)*shp(3,iplot))*shz(iplot) ) * pas
!
        dy(iplot) = ((v(i1,iet  )*shp(1,iplot)
     &              + v(i2,iet  )*shp(2,iplot)
     &              + v(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &              +(v(i1,iet+1)*shp(1,iplot)
     &              + v(i2,iet+1)*shp(2,iplot)
     &              + v(i3,iet+1)*shp(3,iplot))*shz(iplot) ) * pas
!
        IF(sigma) THEN
          deltaz =  (z(i1,iet+1)-z(i1,iet))*shp(1,iplot)
     &            + (z(i2,iet+1)-z(i2,iet))*shp(2,iplot)
     &            + (z(i3,iet+1)-z(i3,iet))*shp(3,iplot)
!
          IF(deltaz.GT.epsdz) THEN
!           DIVISION BY DELTAZ IS DUE TO THE FACT THAT W IS
!           W* MULTIPLIED BY DELTAZ (IT STEMS FROM TRIDW2 IN TELEMAC3D)
            dz(iplot) = ((w(i1,iet  )*shp(1,iplot)
     &                  + w(i2,iet  )*shp(2,iplot)
     &                  + w(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &                  +(w(i1,iet+1)*shp(1,iplot)
     &                  + w(i2,iet+1)*shp(2,iplot)
     &                  + w(i3,iet+1)*shp(3,iplot))*shz(iplot) )
     &                  * pas * (zstar(iet+1)-zstar(iet)) / deltaz
          ELSE
            dz(iplot) = 0.d0
          ENDIF
        ELSE
          dz(iplot) = ((w(i1,iet  )*shp(1,iplot)
     &                + w(i2,iet  )*shp(2,iplot)
     &                + w(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &                +(w(i1,iet+1)*shp(1,iplot)
     &                + w(i2,iet+1)*shp(2,iplot)
     &                + w(i3,iet+1)*shp(3,iplot))*shz(iplot) ) * pas
        ENDIF
!
        xp = xplot(iplot) + dx(iplot)
        yp = yplot(iplot) + dy(iplot)
        zp = zplot(iplot) + dz(iplot)
!
        shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                 -(y(i3)-y(i2))*(xp-x(i2))) * surdet(iel)
        shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                 -(y(i1)-y(i3))*(xp-x(i3))) * surdet(iel)
        shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                 -(y(i2)-y(i1))*(xp-x(i1))) * surdet(iel)
        IF(sigma) THEN
          shz(iplot) = (zp-zstar(iet)) / (zstar(iet+1)-zstar(iet))
        ELSE
          zdown=shp(1,iplot)*z(i1,iet)
     &         +shp(2,iplot)*z(i2,iet)
     &         +shp(3,iplot)*z(i3,iet)
          zup  =shp(1,iplot)*z(i1,iet+1)
     &         +shp(2,iplot)*z(i2,iet+1)
     &         +shp(3,iplot)*z(i3,iet+1)
          shz(iplot) = (zp-zdown) / max(zup-zdown,epsdz)
        ENDIF
!
        xplot(iplot) = xp
        yplot(iplot) = yp
        zplot(iplot) = zp
!
        IF(add) THEN
!         CONTINUOUS SETTING OF THE REACHED POSITION FOR IPLOT
!         AND THE NUMBER OF STEPS DONE ALREADY
          recvchar(iplot)%XP=xplot(iplot)
          recvchar(iplot)%YP=yplot(iplot)
          recvchar(iplot)%ZP=zplot(iplot)
          recvchar(iplot)%DX=dx(iplot)
          recvchar(iplot)%DY=dy(iplot)
          recvchar(iplot)%DZ=dz(iplot)
          recvchar(iplot)%INE=elt(iplot)
          recvchar(iplot)%ISP=isp
        ENDIF
!
!-----------------------------------------------------------------------
!       TEST: IS THE PATHLINE WENT OUT THE ORIGINAL ELEMENT
!-----------------------------------------------------------------------
!
50      CONTINUE
!
        iso = 0
        IF(shp(1,iplot).LT.epsilo)      iso=ibset(iso,2)
        IF(shp(2,iplot).LT.epsilo)      iso=ibset(iso,3)
        IF(shp(3,iplot).LT.epsilo)      iso=ibset(iso,4)
        IF(shz(iplot)  .LT.epsilo)      iso=ibset(iso,0)
        IF(shz(iplot)  .GT.1.d0-epsilo) iso=ibset(iso,1)
!
        IF(iso.NE.0) THEN
!
!-----------------------------------------------------------------------
!         HERE WE ARE OUT OF THE ELEMENT
!-----------------------------------------------------------------------
!
          isoh = iand(iso,28)
          isov = iand(iso, 3)
          iel = elt(iplot)
          iet = eta(iplot)
          xp = xplot(iplot)
          yp = yplot(iplot)
          zp = zplot(iplot)
!
          IF(isoh.NE.0) THEN
!
            IF(isoh.EQ.4) THEN
              ifa = 2
            ELSEIF(isoh.EQ.8) THEN
              ifa = 3
            ELSEIF(isoh.EQ.16) THEN
              ifa = 1
            ELSEIF(isoh.EQ.12) THEN
              ifa = 2
              IF(dx(iplot)*(y(ikle2(iel,3))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,3))-xp)) ifa = 3
            ELSEIF(isoh.EQ.24) THEN
              ifa = 3
              IF(dx(iplot)*(y(ikle2(iel,1))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,1))-xp)) ifa = 1
            ELSE
              ifa = 1
              IF(dx(iplot)*(y(ikle2(iel,2))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,2))-xp)) ifa = 2
            ENDIF
!
            IF(isov.GT.0) THEN
              IF(sigma) THEN
                IF(abs(dz(iplot)).GT.epsdz) THEN
!                 PERCENTAGE OF DISPLACEMENT DONE OUT OF THE ELEMENT
                  a1 = (zp-zstar(iet+isov-1)) / dz(iplot)
                ELSE
                  a1 = 0.d0
                ENDIF
                i1 = ikle2(iel,ifa)
                i2 = ikle2(iel,isui(ifa))
!               IF EXIT POINT THROUGH LEVEL STILL IN TRIANGLE
!               THEN THE REAL EXIT WAS FACES 4 OR 5
!               UPPER AND LOWER TRIANGLE
                IF ((x(i2)-x(i1))*(yp-a1*dy(iplot)-y(i1)).GT.
     &              (y(i2)-y(i1))*(xp-a1*dx(iplot)-x(i1))) ifa=isov+3
              ELSE
                denom=-(x(i2)-x(i1))*dy(iplot)+(y(i2)-y(i1))*dx(iplot)
!               PERCENTAGE OF DISPLACEMENT DONE IN THE ELEMENT
                IF(abs(denom).GT.1.d-10) THEN
                  num=-(xp-x(i1))*dy(iplot)+(yp-y(i1))*dx(iplot)
                  a1=num/denom
                ELSE
                  a1=0.d0
                ENDIF
                zdown=      a1 *z(i2,iet)
     &               +(1.d0-a1)*z(i1,iet)
                zup  =      a1 *z(i2,iet+1)
     &               +(1.d0-a1)*z(i1,iet+1)
!               ZZ: ELEVATION WHEN CROSSING SEGMENT I1-I2
                zz   = zp-(1.d0-a1)*dz(iplot)
!               EXIT THROUGH LOWER OR UPPER TRIANGLE
                IF(zz.GT.zup.OR.zz.LT.zdown) ifa=isov+3
              ENDIF
            ENDIF
!
          ELSE
!
            ifa = isov + 3
!
          ENDIF
!
          iel = ibor(iel,ifa,iet)
!
          IF(ifa.LE.3) THEN
!
!-----------------------------------------------------------------------
!           HERE WE ARRIVE IN ANOTHER ELEMENT THROUGH A QUADRANGULAR FACE
!-----------------------------------------------------------------------
!
            IF(iel.GT.0) THEN
!
!-----------------------------------------------------------------------
!             RELOCALISING IN ADJACENT ELEMENT
!-----------------------------------------------------------------------
!
              i1 = ikle2(iel,1)
              i2 = ikle2(iel,2)
              i3 = ikle2(iel,3)
              elt(iplot) = iel
              shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
              shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
              shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
!
              GOTO 50
!
            ENDIF
!
!-----------------------------------------------------------------------
!           HERE WE PASS TO A NEIGHBOUR SUBDOMAIN AND COLLECT DATA
!-----------------------------------------------------------------------
!
          IF(iel.EQ.-2) THEN
            IF(.NOT.add) THEN
!             INTERFACE CROSSING
              CALL collect_char
     &            (ipid,iplot,elt(iplot),ifa,eta(iplot),0,isp,
     &             nsp,xplot(iplot),yplot(iplot),
     &             zplot(iplot),0.d0,
     &             dx(iplot),dy(iplot),dz(iplot),0.d0,
     &             ifapar,nchdim,nchara)
            ELSE
!             A LOST-AGAIN TRACEBACK DETECTED
!             PROCESSOR NUMBER
              recvchar(iplot)%NEPID=ifapar(ifa,elt(iplot))
              recvchar(iplot)%INE=ifapar(ifa+3,elt(iplot))
              recvchar(iplot)%KNE=eta(iplot)
            ENDIF
!           EXITING LOOP ON ISP
            EXIT
          ENDIF
!
!-----------------------------------------------------------------------
!           SPECIAL TREATMENT FOR SOLID OR LIQUID BOUNDARIES
!-----------------------------------------------------------------------
!
            dxp = dx(iplot)
            dyp = dy(iplot)
            i1  = ikle2(elt(iplot),ifa)
            i2  = ikle2(elt(iplot),isui(ifa))
            dx1 = x(i2) - x(i1)
            dy1 = y(i2) - y(i1)
!
            IF(iel.EQ.-1) THEN
!
!-----------------------------------------------------------------------
!             HERE SOLID BOUNDARY, VELOCITY IS PROJECTED ON THE BOUNDARY
!             AND WE GO ON
!-----------------------------------------------------------------------
!
              a1 = (dxp*dx1+dyp*dy1) / (dx1**2+dy1**2)
              dx(iplot) = a1 * dx1
              dy(iplot) = a1 * dy1
!
              a1=((xp-x(i1))*dx1+(yp-y(i1))*dy1)/(dx1**2+dy1**2)
              shp(          ifa  ,iplot) = 1.d0 - a1
              shp(     isui(ifa) ,iplot) = a1
              shp(isui(isui(ifa)),iplot) = 0.d0
              xplot(iplot) = x(i1) + a1 * dx1
              yplot(iplot) = y(i1) + a1 * dy1
              IF(add) THEN
                recvchar(iplot)%XP=xplot(iplot)
                recvchar(iplot)%YP=yplot(iplot)
                recvchar(iplot)%ZP=zplot(iplot)
                recvchar(iplot)%DX=dx(iplot)
                recvchar(iplot)%DY=dy(iplot)
                recvchar(iplot)%DZ=dz(iplot)
              ENDIF
!
              GOTO 50
!
            ELSEIF(iel.EQ.0) THEN
!
!-----------------------------------------------------------------------
!           HERE WE HAVE A LIQUID BOUNDARY, THE CHARACTERISTIC IS STOPPED
!-----------------------------------------------------------------------
!
            denom = dxp*dy1-dyp*dx1
            IF(abs(denom).GT.1.d-12) THEN
              a1 = (dxp*(yp-y(i1))-dyp*(xp-x(i1))) / denom
            ELSE
              a1 = 0.d0
            ENDIF
            IF(a1.GT.1.d0) a1 = 1.d0
            IF(a1.LT.0.d0) a1 = 0.d0
            shp(          ifa  ,iplot) = 1.d0 - a1
            shp(     isui(ifa) ,iplot) = a1
            shp(isui(isui(ifa)),iplot) = 0.d0
            xplot(iplot) = x(i1) + a1 * dx1
            yplot(iplot) = y(i1) + a1 * dy1
            IF(abs(dxp).GT.abs(dyp)) THEN
              a1 = (xp-xplot(iplot))/dxp
            ELSE
              a1 = (yp-yplot(iplot))/dyp
            ENDIF
            zplot(iplot) = zp - a1*dz(iplot)
            IF(sigma) THEN
              shz(iplot) = (zplot(iplot)-zstar(iet))
     &                   / (zstar(iet+1)-zstar(iet))
            ELSE
              shz(iplot) = (zp-zdown) / max(zup-zdown,epsdz)
            ENDIF
!           THIS IS A MARKER FOR PARTICLES EXITING A DOMAIN
!           SENS=-1 FOR BACKWARD CHARACTERISTICS
            elt(iplot) = - sens * elt(iplot)
!           EXITING LOOP ON ISP
            EXIT
!
            ELSE
!
              WRITE(lu,*) 'UNEXPECTED CASE IN SCHAR41'
              WRITE(lu,*) 'IEL=',iel
              CALL plante(1)
              stop
!
            ENDIF
!
          ELSE
!
!-----------------------------------------------------------------------
!  CASE IFA = 4 OR 5
!  HERE WE EXIT THROUGH TOP OR BOTTOM OF THE PRISM
!-----------------------------------------------------------------------
!
            ifa = ifa - 4
!           HENCE IFA NOW EQUALS 0 OR 1
!
            IF(iel.EQ.1) THEN
!
!-----------------------------------------------------------------------
!           NO NEED TO RECOMPUTE THE VELOCITIES,
!           RELOCALISING IN NEW ELEMENT
!-----------------------------------------------------------------------
!
              eta(iplot) = iet + ifa + ifa - 1
              IF(sigma) THEN
                shz(iplot) = (zp-zstar(eta(iplot)))
     &                     / (zstar(eta(iplot)+1)-zstar(eta(iplot)))
              ELSE
                zdown=shp(1,iplot)*z(i1,eta(iplot))
     &               +shp(2,iplot)*z(i2,eta(iplot))
     &               +shp(3,iplot)*z(i3,eta(iplot))
                zup  =shp(1,iplot)*z(i1,eta(iplot)+1)
     &               +shp(2,iplot)*z(i2,eta(iplot)+1)
     &               +shp(3,iplot)*z(i3,eta(iplot)+1)
                shz(iplot) = (zp-zdown) / max(zup-zdown,epsdz)
              ENDIF
!
              IF(add) THEN
                recvchar(iplot)%KNE=eta(iplot)
              ENDIF
!
              iso = isoh
!
              IF(shz(iplot).LT.     epsilo) iso=ibset(iso,0)
              IF(shz(iplot).GT.1.d0-epsilo) iso=ibset(iso,1)
!
              GOTO 50
!
            ENDIF
!
            IF(iel.EQ.-1) THEN
!
!-----------------------------------------------------------------------
!  LA, ON SAIT QUE LA FACE DE SORTIE EST UNE FRONTIERE SOLIDE
!  ON PROJETTE LE RELIQUAT SUR LA FRONTIERE PUIS ON SE RELOCALISE
!-----------------------------------------------------------------------
!
              dz(iplot) = 0.d0
!
              IF(sigma) THEN
                zplot(iplot) = zstar(iet+ifa)
              ELSE
                zplot(iplot) = shp(1,iplot)*z(i1,iet+ifa)
     &                        +shp(2,iplot)*z(i2,iet+ifa)
     &                        +shp(3,iplot)*z(i3,iet+ifa)
              ENDIF
              shz(iplot) = ifa
!
              iso = isoh
              IF(isoh.NE.0) GOTO 50
!
            ELSE
!
!-----------------------------------------------------------------------
!  LA, ON SAIT QUE LA FACE DE SORTIE EST UNE FRONTIERE LIQUIDE (CAS 0)
!      ON ARRETE ALORS LA REMONTEE DES CARACTERISTIQUES (SIGNE DE ELT)
!  OU, QUE L'ON VIENT DE TRAVERSER UN PLAN AVEC RECALCUL DES VITESSES
!  DEMANDE (CAS 2)
!-----------------------------------------------------------------------
!
              IF(sigma) THEN
                IF(abs(dz(iplot)).GT.epsdz) THEN
                  a1 = (zp-zstar(iet+ifa)) / dz(iplot)
                ELSE
                  a1 = 0.d0
                ENDIF
                xp = xp - a1*dx(iplot)
                yp = yp - a1*dy(iplot)
                zp = zstar(iet+ifa)
              ELSE
                WRITE(lu,*) 'SORTIE EN VERTICALE PAR FRONTIERE LIQUIDE'
                WRITE(lu,*) 'CAS NON PROGRAMME'
                CALL plante(1)
                stop
              ENDIF
              iele = elt(iplot)
              i1 = ikle2(iele,1)
              i2 = ikle2(iele,2)
              i3 = ikle2(iele,3)
!
              shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iele)
              shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iele)
              shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iele)
!
              IF(iel.EQ.2) THEN
!
!-----------------------------------------------------------------------
!  LA, ON SAIT QUE LA FACE DE SORTIE SE SITUE SUR LE PLAN OU ON DEMANDE
!  UN RECALCUL DES VITESSES
!-----------------------------------------------------------------------
!
!               IF IFA = 1 EXIT THROUGH THE TOP
!               IF IFA = 0 EXIT THROUGH THE BOTTOM
!               THEN NEW IET IS  IET+1 IF IFA = 1
!                            AND IET-1 IF IFA = 0
!               THIS IS SUMMARISED BY IET=IET+2*IFA-1
!
!               RECOMPUTED VELOCITIES MUST BE TAKEN AT IET2=IET+IFA
!               I.E. BOTTOM IF EXIT THROUGH THE BOTTOM
!               AND TOP IF EXIT THROUGH THE TOP
!
                iet2 = iet + ifa
                iet  = iet + ifa + ifa - 1
!
                pas2 = pas2 * a1
!
                dx(iplot) = ( u(i1,iet2)*shp(1,iplot)
     &                      + u(i2,iet2)*shp(2,iplot)
     &                      + u(i3,iet2)*shp(3,iplot) ) * pas2
!
                dy(iplot) = ( v(i1,iet2)*shp(1,iplot)
     &                      + v(i2,iet2)*shp(2,iplot)
     &                      + v(i3,iet2)*shp(3,iplot) ) * pas2
!
                IF(sigma) THEN
                  deltaz =  (z(i1,iet+1)-z(i1,iet))*shp(1,iplot)
     &                    + (z(i2,iet+1)-z(i2,iet))*shp(2,iplot)
     &                    + (z(i3,iet+1)-z(i3,iet))*shp(3,iplot)
!
                  IF(deltaz.GT.epsdz) THEN
                    dz(iplot) = ( w(i1,iet2)*shp(1,iplot)
     &                          + w(i2,iet2)*shp(2,iplot)
     &                          + w(i3,iet2)*shp(3,iplot) ) * pas2
     &                          * (zstar(iet+1)-zstar(iet)) / deltaz
                  ELSE
                    dz(iplot) = 0.d0
                  ENDIF
                ELSE
                  dz(iplot)=((w(i1,iet  )*shp(1,iplot)
     &                    + w(i2,iet  )*shp(2,iplot)
     &                    + w(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &                    +(w(i1,iet+1)*shp(1,iplot)
     &                    + w(i2,iet+1)*shp(2,iplot)
     &                    + w(i3,iet+1)*shp(3,iplot))*shz(iplot) )*pas
                ENDIF
!
                xp = xp + dx(iplot)
                yp = yp + dy(iplot)
                zp = zp + dz(iplot)
!
                shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2))) * surdet(iele)
                shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3))) * surdet(iele)
                shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1))) * surdet(iele)
                IF(sigma) THEN
                  shz(iplot)=(zp-zstar(iet))/(zstar(iet+1)-zstar(iet))
                ELSE
                  zdown=shp(1,iplot)*z(i1,iet)
     &                 +shp(2,iplot)*z(i2,iet)
     &                 +shp(3,iplot)*z(i3,iet)
                  zup  =shp(1,iplot)*z(i1,iet+1)
     &                 +shp(2,iplot)*z(i2,iet+1)
     &                 +shp(3,iplot)*z(i3,iet+1)
                  shz(iplot) = (zp-zdown) / max(zup-zdown,epsdz)
                ENDIF
!
                xplot(iplot) = xp
                yplot(iplot) = yp
                zplot(iplot) = zp
                eta(iplot) = iet
!
                iso = 0
!
                IF(shp(1,iplot).LT.epsilo) iso=ibset(iso,2)
                IF(shp(2,iplot).LT.epsilo) iso=ibset(iso,3)
                IF(shp(3,iplot).LT.epsilo) iso=ibset(iso,4)
!
                IF(shz(iplot).LT.     epsilo) iso=ibset(iso,0)
                IF(shz(iplot).GT.1.d0-epsilo) iso=ibset(iso,1)
!
                IF(add) THEN
                  recvchar(iplot)%XP=xplot(iplot)
                  recvchar(iplot)%YP=yplot(iplot)
                  recvchar(iplot)%ZP=zplot(iplot)
                  recvchar(iplot)%DX=dx(iplot)
                  recvchar(iplot)%DY=dy(iplot)
                  recvchar(iplot)%DZ=dz(iplot)
                ENDIF
!
                GOTO 50
!
              ENDIF
!
              xplot(iplot) = xp
              yplot(iplot) = yp
              zplot(iplot) = zp
              shz(iplot) = ifa
!             THIS IS A MARKER FOR PARTICLES EXITING A DOMAIN
!             SENS=-1 FOR BACKWARD CHARACTERISTICS
              elt(iplot) = - sens * elt(iplot)
!             EXITING LOOP ON ISP
              EXIT
!
            ENDIF
!
          ENDIF
!
!       IF(ISO.NE.0)
        ENDIF
!
        ENDDO
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE schar41
!                       **********************
                        SUBROUTINE schar41_sto
!                       **********************
!
     &( u , v , w , dt , nrk , x , y , zstar , z , ikle2 , ibor ,
     &  xplot , yplot , zplot , dx , dy , dz , shp , shz , elt , eta ,
     &  nplot , npoin2 , nelem2 , nelmax2,nplan , surdet ,
     &  sens  , ifapar, nchdim,nchara,add,sigma,viscvi,stocha)
!
!***********************************************************************
! BIEF VERSION 7.0
!
!***********************************************************************
!
!  FONCTION : This is a mere copy of SCHAR41, but with a stochastic
!             model on the horizontal.
!
!-----------------------------------------------------------------------
!                             ARGUMENTS
! .________________.____.______________________________________________.
! |      NOM       |MODE|                   ROLE                       |
! |________________|____|______________________________________________|
! |    U,V,W       | -->| COMPONENTS OF ADVECTION VELOCITY
! |                |    | BUT W IS W* x DELTAZ (STEMS FROM TRIDW2)
! |    DT          | -->| PAS DE TEMPS.                                |
! |    NRK         | -->| NOMBRE DE SOUS-PAS DE RUNGE-KUTTA.           |
! |    X,Y,ZSTAR   | -->| COORDONNEES DES POINTS DU MAILLAGE.          |
! |    Z           | -->| COTE DANS LE MAILLAGE REEL                   |
! |    IKLE2       | -->| TRANSITION ENTRE LES NUMEROTATIONS LOCALE    |
! |                |    | ET GLOBALE DU MAILLAGE 2D.                   |
! |    IBOR        | -->| NUMEROS 2D DES ELEMENTS AYANT UNE FACE COMMUNE
! |                |    | AVEC L'ELEMENT .  SI IFABOR<0 OU NUL         |
! |                |    | ON A UNE FACE LIQUIDE,SOLIDE,OU PERIODIQUE   |
! |  X..,Y..,ZPLOT |<-->| POSITIONS SUCCESSIVES DES DERIVANTS.         |
! |    DX,DY,DZ    | -- | STOCKAGE DES SOUS-PAS .                      |
! |    SHP         |<-->| COORDONNEES BARYCENTRIQUES 2D AU PIED DES    |
! |                |    | COURBES CARACTERISTIQUES.                    |
! |    SHZ         |<-->| COORDONNEES BARYCENTRIQUES SUIVANT Z DES     |
! |                |    | NOEUDS DANS LEURS ETAGES "ETA" ASSOCIES.     |
! |    ELT         |<-->| NUMEROS DES ELEMENTS 2D CHOISIS POUR CHAQUE  |
! |                |    | NOEUD.                                       |
! |    ETA         |<-->| NUMEROS DES ETAGES CHOISIS POUR CHAQUE NOEUD.|
! |    NPLOT       | -->| NOMBRE DE DERIVANTS.                         |
! |    NPOIN2      | -->| NOMBRE DE POINTS DU MAILLAGE 2D.             |
! |    NELEM2      | -->| NOMBRE D'ELEMENTS DU MAILLAGE 2D.            |
! |    NPLAN       | -->| NOMBRE DE PLANS.
! |    SIGMA       | -->| IF YES, TRANSFORMED MESH
! |    SURDET      | -->| VARIABLE UTILISEE PAR LA TRANSFORMEE ISOPARAM.
! |    SENS        | -->| DESCENTE OU REMONTEE DES CARACTERISTIQUES.   |
! |    ISO         | -- | STOCKAGE BINAIRE DE LA FACE DE SORTIE.       |
! |________________|____|______________________________________________|
!  MODE: -->(DONNEE NON MODIFIEE),<--(RESULTAT),<-->(DONNEE MODIFIEE)
!-----------------------------------------------------------------------
!     - APPELE PAR : CARACT
!     - PROGRAMMES APPELES : NEANT
!
!***********************************************************************
!
      USE bief
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)    :: SENS,NPLAN,NCHDIM,STOCHA
      INTEGER         , INTENT(IN)    :: NPOIN2,NELEM2,NPLOT,NRK,NELMAX2
      INTEGER         , INTENT(IN)    :: IKLE2(nelmax2,3)
      INTEGER         , INTENT(INOUT) :: ELT(nplot),NCHARA
      DOUBLE PRECISION, INTENT(IN)    :: U(npoin2,nplan),V(npoin2,nplan)
      DOUBLE PRECISION, INTENT(IN)    :: W(npoin2,nplan),SURDET(nelem2)
      DOUBLE PRECISION, INTENT(INOUT) :: XPLOT(nplot),YPLOT(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: ZPLOT(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: SHP(3,nplot),SHZ(nplot)
      DOUBLE PRECISION, INTENT(IN)    :: X(npoin2),Y(npoin2),DT
      DOUBLE PRECISION, INTENT(IN)    :: Z(npoin2,nplan),ZSTAR(nplan)
      DOUBLE PRECISION, INTENT(INOUT) :: DX(nplot),DY(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: DZ(nplot)
      INTEGER         , INTENT(IN)    :: IBOR(nelmax2,5,nplan-1)
      INTEGER         , INTENT(INOUT) :: ETA(nplot)
      INTEGER         , INTENT(IN)    :: IFAPAR(6,*)
      LOGICAL, INTENT(IN)             :: ADD,SIGMA
      type(bief_obj), INTENT(INOUT)  :: viscvi
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELE,ISO,ISPDONE,NSP
      INTEGER :: IPLOT,ISP,I1,I2,I3,IEL,IET,IET2,ISOH,ISOV,IFA,ISUI(3)
!
      DOUBLE PRECISION PAS,A1,DX1,DY1,DXP,DYP,XP,YP,ZP,NUM,DENOM
      DOUBLE PRECISION DELTAZ,PAS2,ZUP,ZDOWN,ZZ
!
      INTRINSIC abs , int , max , sqrt
!
!     FOR STOCHASTIC DIFFUSION
      DOUBLE PRECISION RAND1,RAND2,A,B,C,D,E,DIFF_X,DIFF_Y,DEUXPI
!
      parameter( isui = (/ 2 , 3 , 1 /) )
      DOUBLE PRECISION, PARAMETER :: EPSILO = -1.d-6
      DOUBLE PRECISION, PARAMETER :: EPSDZ  = 1.d-4
!
!-----------------------------------------------------------------------
!     FOR EVERY POINT
!-----------------------------------------------------------------------
!
      DO iplot = 1 , nplot
!
        IF(add) THEN
!
        xplot(iplot)   = recvchar(iplot)%XP
        yplot(iplot)   = recvchar(iplot)%YP
        zplot(iplot)   = recvchar(iplot)%ZP
        dx(iplot)      = recvchar(iplot)%DX
        dy(iplot)      = recvchar(iplot)%DY
        dz(iplot)      = recvchar(iplot)%DZ
        elt(iplot)     = recvchar(iplot)%INE
        eta(iplot)     = recvchar(iplot)%KNE
        nsp            = recvchar(iplot)%NSP ! R-K STEPS TO BE FULLFILLED
        ispdone        = recvchar(iplot)%ISP ! R-K STEPS ALREADY DONE
        iel = elt(iplot)
        iet = eta(iplot)
        xp  = xplot(iplot)
        yp  = yplot(iplot)
        zp  = zplot(iplot)
        i1 = ikle2(iel,1)
        i2 = ikle2(iel,2)
        i3 = ikle2(iel,3)
        shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                 -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
        shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                 -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
        shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                 -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
        IF(sigma) THEN
          shz(iplot) = (zp-zstar(iet)) / (zstar(iet+1)-zstar(iet))
        ELSE
          zdown=shp(1,iplot)*z(i1,iet)
     &         +shp(2,iplot)*z(i2,iet)
     &         +shp(3,iplot)*z(i3,iet)
          zup  =shp(1,iplot)*z(i1,iet+1)
     &         +shp(2,iplot)*z(i2,iet+1)
     &         +shp(3,iplot)*z(i3,iet+1)
          shz(iplot) = (zp-zdown) / max(zup-zdown,epsdz)
        ENDIF
!       ASSUME ALL ARE LOCALISED, IT WILL BE SET OTHERWISE IF LOST-AGAIN
        recvchar(iplot)%NEPID=-1
!
        ELSE
!
          iel = elt(iplot)
!         POINTS WITH IEL=0 ARE TREATED SO THAT THE FINAL
!         INTERPOLATION GIVES 0.,
!         AND WE SKIP TO NEXT POINT IPLOT (CYCLE)
          IF(iel.EQ.0) THEN
            elt(iplot)=1
            eta(iplot)=1
            shp(1,iplot)=0.d0
            shp(2,iplot)=0.d0
            shp(3,iplot)=0.d0
            shz(iplot)=0.d0
            cycle
          ENDIF
          iet = eta(iplot)
          i1 = ikle2(iel,1)
          i2 = ikle2(iel,2)
          i3 = ikle2(iel,3)
          dxp =( u(i1,iet  )*shp(1,iplot)
     &         + u(i2,iet  )*shp(2,iplot)
     &         + u(i3,iet  )*shp(3,iplot) )*(1.d0-shz(iplot))
     &        +( u(i1,iet+1)*shp(1,iplot)
     &         + u(i2,iet+1)*shp(2,iplot)
     &         + u(i3,iet+1)*shp(3,iplot) )*shz(iplot)
          dyp =( v(i1,iet  )*shp(1,iplot)
     &         + v(i2,iet  )*shp(2,iplot)
     &         + v(i3,iet  )*shp(3,iplot) )*(1.d0-shz(iplot))
     &        +( v(i1,iet+1)*shp(1,iplot)
     &         + v(i2,iet+1)*shp(2,iplot)
     &         + v(i3,iet+1)*shp(3,iplot) )*shz(iplot)
!         VERTICAL VELOCITY NOT CONSIDERED HERE !!
          nsp=max(1,int(nrk*dt*sqrt((dxp**2+dyp**2)*surdet(iel))))
          ispdone=1
!
        ENDIF
!
        pas = sens * dt / nsp
!
!       LOOP ON RUNGE-KUTTA SUB-STEPS
!
!       COMPILER MUST DO NOTHING IF ISPDONE>NSP
!       IN MODE "ADD", ISP = ISPDONE HAS NOT BEEN FULLY DONE
!       IT IS RESTARTED HERE
!
      DO isp = ispdone,nsp
!
!-----------------------------------------------------------------------
!       LOCALISING THE ARRIVAL POINT
!-----------------------------------------------------------------------
!
        pas2 = pas
!
!       IN MODE "ADD" ITERATIONS ALREADY DONE ARE SKIPPED AND
!                     CHARACTERISTICS GONE IN ANOTHER SUB-DOMAIN SKIPPED
!
        IF(add) THEN
          IF(isp.EQ.ispdone) GO TO 50
          IF(recvchar(iplot)%NEPID.NE.-1) cycle
        ENDIF
!
        iel = elt(iplot)
        iet = eta(iplot)
        i1 = ikle2(iel,1)
        i2 = ikle2(iel,2)
        i3 = ikle2(iel,3)
!
        dx(iplot) = ((u(i1,iet  )*shp(1,iplot)
     &              + u(i2,iet  )*shp(2,iplot)
     &              + u(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &              +(u(i1,iet+1)*shp(1,iplot)
     &              + u(i2,iet+1)*shp(2,iplot)
     &              + u(i3,iet+1)*shp(3,iplot))*shz(iplot) ) * pas
!
        dy(iplot) = ((v(i1,iet  )*shp(1,iplot)
     &              + v(i2,iet  )*shp(2,iplot)
     &              + v(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &              +(v(i1,iet+1)*shp(1,iplot)
     &              + v(i2,iet+1)*shp(2,iplot)
     &              + v(i3,iet+1)*shp(3,iplot))*shz(iplot) ) * pas
!
        IF(sigma) THEN
          deltaz =  (z(i1,iet+1)-z(i1,iet))*shp(1,iplot)
     &            + (z(i2,iet+1)-z(i2,iet))*shp(2,iplot)
     &            + (z(i3,iet+1)-z(i3,iet))*shp(3,iplot)
!
          IF(deltaz.GT.epsdz) THEN
!           DIVISION BY DELTAZ IS DUE TO THE FACT THAT W IS
!           W* MULTIPLIED BY DELTAZ (IT STEMS FROM TRIDW2 IN TELEMAC3D)
            dz(iplot) = ((w(i1,iet  )*shp(1,iplot)
     &                  + w(i2,iet  )*shp(2,iplot)
     &                  + w(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &                  +(w(i1,iet+1)*shp(1,iplot)
     &                  + w(i2,iet+1)*shp(2,iplot)
     &                  + w(i3,iet+1)*shp(3,iplot))*shz(iplot) )
     &                  * pas * (zstar(iet+1)-zstar(iet)) / deltaz
          ELSE
            dz(iplot) = 0.d0
          ENDIF
        ELSE
          dz(iplot) = ((w(i1,iet  )*shp(1,iplot)
     &                + w(i2,iet  )*shp(2,iplot)
     &                + w(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &                +(w(i1,iet+1)*shp(1,iplot)
     &                + w(i2,iet+1)*shp(2,iplot)
     &                + w(i3,iet+1)*shp(3,iplot))*shz(iplot) ) * pas
        ENDIF
!
!       STOCHASTIC DIFFUSION
!
        IF(stocha.EQ.1) THEN
!         COMPUTING LOCAL VISCOSITY
          a=max(viscvi%ADR(1)%P%R(i1)*shp(1,iplot)
     &         +viscvi%ADR(1)%P%R(i2)*shp(2,iplot)
     &         +viscvi%ADR(1)%P%R(i3)*shp(3,iplot),0.d0)
          b=max(viscvi%ADR(2)%P%R(i1)*shp(1,iplot)
     &         +viscvi%ADR(2)%P%R(i2)*shp(2,iplot)
     &         +viscvi%ADR(2)%P%R(i3)*shp(3,iplot),0.d0)
!         DISPLACEMENT DUE TO RANDOM DIFFUSION
          deuxpi=2.d0*acos(-1.d0)
          CALL random_number(rand1)
          CALL random_number(rand2)
          c=sqrt(-2.d0*log(rand1))
          d=c*cos(deuxpi*rand2)
          e=c*sin(deuxpi*rand2)
          diff_x=d*sqrt(2.d0*a/0.72d0)
          diff_y=e*sqrt(2.d0*b/0.72d0)
          dx(iplot) = dx(iplot) + diff_x*sqrt(abs(pas))
          dy(iplot) = dy(iplot) + diff_y*sqrt(abs(pas))
        ENDIF

        xp = xplot(iplot) + dx(iplot)
        yp = yplot(iplot) + dy(iplot)
        zp = zplot(iplot) + dz(iplot)
!
        shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                 -(y(i3)-y(i2))*(xp-x(i2))) * surdet(iel)
        shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                 -(y(i1)-y(i3))*(xp-x(i3))) * surdet(iel)
        shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                 -(y(i2)-y(i1))*(xp-x(i1))) * surdet(iel)
        IF(sigma) THEN
          shz(iplot) = (zp-zstar(iet)) / (zstar(iet+1)-zstar(iet))
        ELSE
          zdown=shp(1,iplot)*z(i1,iet)
     &         +shp(2,iplot)*z(i2,iet)
     &         +shp(3,iplot)*z(i3,iet)
          zup  =shp(1,iplot)*z(i1,iet+1)
     &         +shp(2,iplot)*z(i2,iet+1)
     &         +shp(3,iplot)*z(i3,iet+1)
          shz(iplot) = (zp-zdown) / max(zup-zdown,epsdz)
        ENDIF
!
        xplot(iplot) = xp
        yplot(iplot) = yp
        zplot(iplot) = zp
!
        IF(add) THEN
!         CONTINUOUS SETTING OF THE REACHED POSITION FOR IPLOT
!         AND THE NUMBER OF STEPS DONE ALREADY
          recvchar(iplot)%XP=xplot(iplot)
          recvchar(iplot)%YP=yplot(iplot)
          recvchar(iplot)%ZP=zplot(iplot)
          recvchar(iplot)%DX=dx(iplot)
          recvchar(iplot)%DY=dy(iplot)
          recvchar(iplot)%DZ=dz(iplot)
          recvchar(iplot)%INE=elt(iplot)
          recvchar(iplot)%ISP=isp
        ENDIF
!
!-----------------------------------------------------------------------
!       TEST: IS THE PATHLINE WENT OUT THE ORIGINAL ELEMENT
!-----------------------------------------------------------------------
!
50      CONTINUE
!
        iso = 0
        IF(shp(1,iplot).LT.epsilo)      iso=ibset(iso,2)
        IF(shp(2,iplot).LT.epsilo)      iso=ibset(iso,3)
        IF(shp(3,iplot).LT.epsilo)      iso=ibset(iso,4)
        IF(shz(iplot)  .LT.epsilo)      iso=ibset(iso,0)
        IF(shz(iplot)  .GT.1.d0-epsilo) iso=ibset(iso,1)
!
        IF(iso.NE.0) THEN
!
!-----------------------------------------------------------------------
!         HERE WE ARE OUT OF THE ELEMENT
!-----------------------------------------------------------------------
!
          isoh = iand(iso,28)
          isov = iand(iso, 3)
          iel = elt(iplot)
          iet = eta(iplot)
          xp = xplot(iplot)
          yp = yplot(iplot)
          zp = zplot(iplot)
!
          IF(isoh.NE.0) THEN
!
            IF(isoh.EQ.4) THEN
              ifa = 2
            ELSEIF(isoh.EQ.8) THEN
              ifa = 3
            ELSEIF(isoh.EQ.16) THEN
              ifa = 1
            ELSEIF(isoh.EQ.12) THEN
              ifa = 2
              IF(dx(iplot)*(y(ikle2(iel,3))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,3))-xp)) ifa = 3
            ELSEIF(isoh.EQ.24) THEN
              ifa = 3
              IF(dx(iplot)*(y(ikle2(iel,1))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,1))-xp)) ifa = 1
            ELSE
              ifa = 1
              IF(dx(iplot)*(y(ikle2(iel,2))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,2))-xp)) ifa = 2
            ENDIF
!
            IF(isov.GT.0) THEN
              IF(sigma) THEN
                IF(abs(dz(iplot)).GT.epsdz) THEN
!                 PERCENTAGE OF DISPLACEMENT DONE OUT OF THE ELEMENT
                  a1 = (zp-zstar(iet+isov-1)) / dz(iplot)
                ELSE
                  a1 = 0.d0
                ENDIF
                i1 = ikle2(iel,ifa)
                i2 = ikle2(iel,isui(ifa))
!               IF EXIT POINT THROUGH LEVEL STILL IN TRIANGLE
!               THEN THE REAL EXIT WAS FACES 4 OR 5
!               UPPER AND LOWER TRIANGLE
                IF ((x(i2)-x(i1))*(yp-a1*dy(iplot)-y(i1)).GT.
     &              (y(i2)-y(i1))*(xp-a1*dx(iplot)-x(i1))) ifa=isov+3
              ELSE
                denom=-(x(i2)-x(i1))*dy(iplot)+(y(i2)-y(i1))*dx(iplot)
!               PERCENTAGE OF DISPLACEMENT DONE IN THE ELEMENT
                IF(abs(denom).GT.1.d-10) THEN
                  num=-(xp-x(i1))*dy(iplot)+(yp-y(i1))*dx(iplot)
                  a1=num/denom
                ELSE
                  a1=0.d0
                ENDIF
                zdown=      a1 *z(i2,iet)
     &               +(1.d0-a1)*z(i1,iet)
                zup  =      a1 *z(i2,iet+1)
     &               +(1.d0-a1)*z(i1,iet+1)
!               ZZ: ELEVATION WHEN CROSSING SEGMENT I1-I2
                zz   = zp-(1.d0-a1)*dz(iplot)
!               EXIT THROUGH LOWER OR UPPER TRIANGLE
                IF(zz.GT.zup.OR.zz.LT.zdown) ifa=isov+3
              ENDIF
            ENDIF
!
          ELSE
!
            ifa = isov + 3
!
          ENDIF
!
          iel = ibor(iel,ifa,iet)
!
          IF(ifa.LE.3) THEN
!
!-----------------------------------------------------------------------
!           HERE WE ARRIVE IN ANOTHER ELEMENT THROUGH A QUADRANGULAR FACE
!-----------------------------------------------------------------------
!
            IF(iel.GT.0) THEN
!
!-----------------------------------------------------------------------
!             RELOCALISING IN ADJACENT ELEMENT
!-----------------------------------------------------------------------
!
              i1 = ikle2(iel,1)
              i2 = ikle2(iel,2)
              i3 = ikle2(iel,3)
              elt(iplot) = iel
              shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
              shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
              shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
!
              GOTO 50
!
            ENDIF
!
!-----------------------------------------------------------------------
!           HERE WE PASS TO A NEIGHBOUR SUBDOMAIN AND COLLECT DATA
!-----------------------------------------------------------------------
!
          IF(iel.EQ.-2) THEN
            IF(.NOT.add) THEN
!             INTERFACE CROSSING
              CALL collect_char
     &            (ipid,iplot,elt(iplot),ifa,eta(iplot),0,isp,
     &             nsp,xplot(iplot),yplot(iplot),
     &             zplot(iplot),0.d0,
     &             dx(iplot),dy(iplot),dz(iplot),0.d0,
     &             ifapar,nchdim,nchara)
            ELSE
!             A LOST-AGAIN TRACEBACK DETECTED
!             PROCESSOR NUMBER
              recvchar(iplot)%NEPID=ifapar(ifa,elt(iplot))
              recvchar(iplot)%INE=ifapar(ifa+3,elt(iplot))
              recvchar(iplot)%KNE=eta(iplot)
            ENDIF
!           EXITING LOOP ON ISP
            EXIT
          ENDIF
!
!-----------------------------------------------------------------------
!           SPECIAL TREATMENT FOR SOLID OR LIQUID BOUNDARIES
!-----------------------------------------------------------------------
!
            dxp = dx(iplot)
            dyp = dy(iplot)
            i1  = ikle2(elt(iplot),ifa)
            i2  = ikle2(elt(iplot),isui(ifa))
            dx1 = x(i2) - x(i1)
            dy1 = y(i2) - y(i1)
!
            IF(iel.EQ.-1) THEN
!
!-----------------------------------------------------------------------
!             HERE SOLID BOUNDARY, VELOCITY IS PROJECTED ON THE BOUNDARY
!             AND WE GO ON
!-----------------------------------------------------------------------
!
              a1 = (dxp*dx1+dyp*dy1) / (dx1**2+dy1**2)
              dx(iplot) = a1 * dx1
              dy(iplot) = a1 * dy1
!
              a1=((xp-x(i1))*dx1+(yp-y(i1))*dy1)/(dx1**2+dy1**2)
              shp(          ifa  ,iplot) = 1.d0 - a1
              shp(     isui(ifa) ,iplot) = a1
              shp(isui(isui(ifa)),iplot) = 0.d0
              xplot(iplot) = x(i1) + a1 * dx1
              yplot(iplot) = y(i1) + a1 * dy1
              IF(add) THEN
                recvchar(iplot)%XP=xplot(iplot)
                recvchar(iplot)%YP=yplot(iplot)
                recvchar(iplot)%ZP=zplot(iplot)
                recvchar(iplot)%DX=dx(iplot)
                recvchar(iplot)%DY=dy(iplot)
                recvchar(iplot)%DZ=dz(iplot)
              ENDIF
!
              GOTO 50
!
            ELSEIF(iel.EQ.0) THEN
!
!-----------------------------------------------------------------------
!           HERE WE HAVE A LIQUID BOUNDARY, THE CHARACTERISTIC IS STOPPED
!-----------------------------------------------------------------------
!
            denom = dxp*dy1-dyp*dx1
            IF(abs(denom).GT.1.d-12) THEN
              a1 = (dxp*(yp-y(i1))-dyp*(xp-x(i1))) / denom
            ELSE
              a1 = 0.d0
            ENDIF
            IF(a1.GT.1.d0) a1 = 1.d0
            IF(a1.LT.0.d0) a1 = 0.d0
            shp(          ifa  ,iplot) = 1.d0 - a1
            shp(     isui(ifa) ,iplot) = a1
            shp(isui(isui(ifa)),iplot) = 0.d0
            xplot(iplot) = x(i1) + a1 * dx1
            yplot(iplot) = y(i1) + a1 * dy1
            IF(abs(dxp).GT.abs(dyp)) THEN
              a1 = (xp-xplot(iplot))/dxp
            ELSE
              a1 = (yp-yplot(iplot))/dyp
            ENDIF
            zplot(iplot) = zp - a1*dz(iplot)
            IF(sigma) THEN
              shz(iplot) = (zplot(iplot)-zstar(iet))
     &                   / (zstar(iet+1)-zstar(iet))
            ELSE
              shz(iplot) = (zp-zdown) / max(zup-zdown,epsdz)
            ENDIF
!           THIS IS A MARKER FOR PARTICLES EXITING A DOMAIN
!           SENS=-1 FOR BACKWARD CHARACTERISTICS
            elt(iplot) = - sens * elt(iplot)
!           EXITING LOOP ON ISP
            EXIT
!
            ELSE
!
              WRITE(lu,*) 'UNEXPECTED CASE IN SCHAR41'
              WRITE(lu,*) 'IEL=',iel
              CALL plante(1)
              stop
!
            ENDIF
!
          ELSE
!
!-----------------------------------------------------------------------
!  CASE IFA = 4 OR 5
!  HERE WE EXIT THROUGH TOP OR BOTTOM OF THE PRISM
!-----------------------------------------------------------------------
!
            ifa = ifa - 4
!           HENCE IFA NOW EQUALS 0 OR 1
!
            IF(iel.EQ.1) THEN
!
!-----------------------------------------------------------------------
!           NO NEED TO RECOMPUTE THE VELOCITIES,
!           RELOCALISING IN NEW ELEMENT
!-----------------------------------------------------------------------
!
              eta(iplot) = iet + ifa + ifa - 1
              IF(sigma) THEN
                shz(iplot) = (zp-zstar(eta(iplot)))
     &                     / (zstar(eta(iplot)+1)-zstar(eta(iplot)))
              ELSE
                zdown=shp(1,iplot)*z(i1,eta(iplot))
     &               +shp(2,iplot)*z(i2,eta(iplot))
     &               +shp(3,iplot)*z(i3,eta(iplot))
                zup  =shp(1,iplot)*z(i1,eta(iplot)+1)
     &               +shp(2,iplot)*z(i2,eta(iplot)+1)
     &               +shp(3,iplot)*z(i3,eta(iplot)+1)
                shz(iplot) = (zp-zdown) / max(zup-zdown,epsdz)
              ENDIF
!
              IF(add) THEN
                recvchar(iplot)%KNE=eta(iplot)
              ENDIF
!
              iso = isoh
!
              IF(shz(iplot).LT.     epsilo) iso=ibset(iso,0)
              IF(shz(iplot).GT.1.d0-epsilo) iso=ibset(iso,1)
!
              GOTO 50
!
            ENDIF
!
            IF(iel.EQ.-1) THEN
!
!-----------------------------------------------------------------------
!  LA, ON SAIT QUE LA FACE DE SORTIE EST UNE FRONTIERE SOLIDE
!  ON PROJETTE LE RELIQUAT SUR LA FRONTIERE PUIS ON SE RELOCALISE
!-----------------------------------------------------------------------
!
              dz(iplot) = 0.d0
!
              IF(sigma) THEN
                zplot(iplot) = zstar(iet+ifa)
              ELSE
                zplot(iplot) = shp(1,iplot)*z(i1,iet+ifa)
     &                        +shp(2,iplot)*z(i2,iet+ifa)
     &                        +shp(3,iplot)*z(i3,iet+ifa)
              ENDIF
              shz(iplot) = ifa
!
              iso = isoh
              IF(isoh.NE.0) GOTO 50
!
            ELSE
!
!-----------------------------------------------------------------------
!  LA, ON SAIT QUE LA FACE DE SORTIE EST UNE FRONTIERE LIQUIDE (CAS 0)
!      ON ARRETE ALORS LA REMONTEE DES CARACTERISTIQUES (SIGNE DE ELT)
!  OU, QUE L'ON VIENT DE TRAVERSER UN PLAN AVEC RECALCUL DES VITESSES
!  DEMANDE (CAS 2)
!-----------------------------------------------------------------------
!
              IF(sigma) THEN
                IF(abs(dz(iplot)).GT.epsdz) THEN
                  a1 = (zp-zstar(iet+ifa)) / dz(iplot)
                ELSE
                  a1 = 0.d0
                ENDIF
                xp = xp - a1*dx(iplot)
                yp = yp - a1*dy(iplot)
                zp = zstar(iet+ifa)
              ELSE
                WRITE(lu,*) 'VERTICAL GETTING OUT FROM LIQUID '
                WRITE(lu,*) 'BOUNDARY - NON IMPLEMENTED CASE'
                CALL plante(1)
                stop
              ENDIF
              iele = elt(iplot)
              i1 = ikle2(iele,1)
              i2 = ikle2(iele,2)
              i3 = ikle2(iele,3)
!
              shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iele)
              shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iele)
              shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iele)
!
              IF(iel.EQ.2) THEN
!
!-----------------------------------------------------------------------
!  LA, ON SAIT QUE LA FACE DE SORTIE SE SITUE SUR LE PLAN OU ON DEMANDE
!  UN RECALCUL DES VITESSES
!-----------------------------------------------------------------------
!
!               IF IFA = 1 EXIT THROUGH THE TOP
!               IF IFA = 0 EXIT THROUGH THE BOTTOM
!               THEN NEW IET IS  IET+1 IF IFA = 1
!                            AND IET-1 IF IFA = 0
!               THIS IS SUMMARISED BY IET=IET+2*IFA-1
!
!               RECOMPUTED VELOCITIES MUST BE TAKEN AT IET2=IET+IFA
!               I.E. BOTTOM IF EXIT THROUGH THE BOTTOM
!               AND TOP IF EXIT THROUGH THE TOP
!
                iet2 = iet + ifa
                iet  = iet + ifa + ifa - 1
!
                pas2 = pas2 * a1
!
                dx(iplot) = ( u(i1,iet2)*shp(1,iplot)
     &                      + u(i2,iet2)*shp(2,iplot)
     &                      + u(i3,iet2)*shp(3,iplot) ) * pas2
!
                dy(iplot) = ( v(i1,iet2)*shp(1,iplot)
     &                      + v(i2,iet2)*shp(2,iplot)
     &                      + v(i3,iet2)*shp(3,iplot) ) * pas2
!
                IF(sigma) THEN
                  deltaz =  (z(i1,iet+1)-z(i1,iet))*shp(1,iplot)
     &                    + (z(i2,iet+1)-z(i2,iet))*shp(2,iplot)
     &                    + (z(i3,iet+1)-z(i3,iet))*shp(3,iplot)
!
                  IF(deltaz.GT.epsdz) THEN
                    dz(iplot) = ( w(i1,iet2)*shp(1,iplot)
     &                          + w(i2,iet2)*shp(2,iplot)
     &                          + w(i3,iet2)*shp(3,iplot) ) * pas2
     &                          * (zstar(iet+1)-zstar(iet)) / deltaz
                  ELSE
                    dz(iplot) = 0.d0
                  ENDIF
                ELSE
                  dz(iplot)=((w(i1,iet  )*shp(1,iplot)
     &                    + w(i2,iet  )*shp(2,iplot)
     &                    + w(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &                    +(w(i1,iet+1)*shp(1,iplot)
     &                    + w(i2,iet+1)*shp(2,iplot)
     &                    + w(i3,iet+1)*shp(3,iplot))*shz(iplot) )*pas
                ENDIF
!
                xp = xp + dx(iplot)
                yp = yp + dy(iplot)
                zp = zp + dz(iplot)
!
                shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2))) * surdet(iele)
                shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3))) * surdet(iele)
                shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1))) * surdet(iele)
                IF(sigma) THEN
                  shz(iplot)=(zp-zstar(iet))/(zstar(iet+1)-zstar(iet))
                ELSE
                  zdown=shp(1,iplot)*z(i1,iet)
     &                 +shp(2,iplot)*z(i2,iet)
     &                 +shp(3,iplot)*z(i3,iet)
                  zup  =shp(1,iplot)*z(i1,iet+1)
     &                 +shp(2,iplot)*z(i2,iet+1)
     &                 +shp(3,iplot)*z(i3,iet+1)
                  shz(iplot) = (zp-zdown) / max(zup-zdown,epsdz)
                ENDIF
!
                xplot(iplot) = xp
                yplot(iplot) = yp
                zplot(iplot) = zp
                eta(iplot) = iet
!
                iso = 0
!
                IF(shp(1,iplot).LT.epsilo) iso=ibset(iso,2)
                IF(shp(2,iplot).LT.epsilo) iso=ibset(iso,3)
                IF(shp(3,iplot).LT.epsilo) iso=ibset(iso,4)
!
                IF(shz(iplot).LT.     epsilo) iso=ibset(iso,0)
                IF(shz(iplot).GT.1.d0-epsilo) iso=ibset(iso,1)
!
                IF(add) THEN
                  recvchar(iplot)%XP=xplot(iplot)
                  recvchar(iplot)%YP=yplot(iplot)
                  recvchar(iplot)%ZP=zplot(iplot)
                  recvchar(iplot)%DX=dx(iplot)
                  recvchar(iplot)%DY=dy(iplot)
                  recvchar(iplot)%DZ=dz(iplot)
                ENDIF
!
                GOTO 50
!
              ENDIF
!
              xplot(iplot) = xp
              yplot(iplot) = yp
              zplot(iplot) = zp
              shz(iplot) = ifa
!             THIS IS A MARKER FOR PARTICLES EXITING A DOMAIN
!             SENS=-1 FOR BACKWARD CHARACTERISTICS
              elt(iplot) = - sens * elt(iplot)
!             EXITING LOOP ON ISP
              EXIT
!
            ENDIF
!
          ENDIF
!
!       IF(ISO.NE.0)
        ENDIF
!
        ENDDO
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE schar41_sto
!                       ************************
                        SUBROUTINE schar41_sigma
!                       ************************
!
     &( u , v , w , dt , nrk , x , y , zstar , z , ikle2 , ibor ,
     &  xplot , yplot , zplot , dx , dy , dz , shp , shz , elt , eta ,
     &  nplot , npoin2 , nelem2 , nelmax2,nplan , surdet ,
     &  sens  , ifapar, nchdim,nchara,add)
!
!***********************************************************************
! BIEF VERSION 7.0           28/04/93     J-M JANIN (LNH) 30 87 72 84
!                        12/10/05     J-M HERVOUET (LNHE) 01 30 87 80 18
!
! 08/11/04 : ADAPTATION A LA TRANSFORMEE SIGMA GENERALISEE
! 12/10/05 : BUG CORRIGE, VOIR VARIABLE IELE QUI ETAIT AVANT IEL
!            ET EFFACAIT UN AUTRE IEL.
!
!
!***********************************************************************
!
!  FONCTION : Exactly like SCHAR41 but optimised for SIGMA=.TRUE.
!             Could be replaced by SCHAR41.
!
!     REMONTEE OU DESCENTE
!     DES COURBES CARACTERISTIQUES
!     SUR DES PRISMES DE TELEMAC-3D
!     DANS L'INTERVALLE DE TEMPS DT
!     AVEC UNE DISCRETISATION ELEMENTS FINIS
!
!
!  DISCRETISATION :
!
!     LE DOMAINE EST APPROCHE PAR UNE DISCRETISATION ELEMENTS FINIS
!     UNE APPROXIMATION LOCALE EST DEFINIE POUR LE VECTEUR VITESSE :
!     LA VALEUR EN UN POINT D'UN ELEMENT NE DEPEND QUE DES VALEURS
!     AUX NOEUDS DE CET ELEMENT
!
!
!  RESTRICTIONS ET HYPOTHESES :
!
!     LE CHAMP CONVECTEUR U EST SUPPOSE INDEPENDANT DU TEMPS
!
!-----------------------------------------------------------------------
!                             ARGUMENTS
! .________________.____.______________________________________________.
! |      NOM       |MODE|                   ROLE                       |
! |________________|____|______________________________________________|
! |    U,V,W       | -->| COMPONENTS OF ADVECTION VELOCITY
! |                |    | BUT W IS W* x DELTAZ (STEMS FROM TRIDW2)
! |    DT          | -->| PAS DE TEMPS.                                |
! |    NRK         | -->| NOMBRE DE SOUS-PAS DE RUNGE-KUTTA.           |
! |    X,Y,ZSTAR   | -->| COORDONNEES DES POINTS DU MAILLAGE.          |
! |    Z           | -->| COTE DANS LE MAILLAGE REEL                   |
! |    IKLE2       | -->| TRANSITION ENTRE LES NUMEROTATIONS LOCALE    |
! |                |    | ET GLOBALE DU MAILLAGE 2D.                   |
! |    IBOR        | -->| NUMEROS 2D DES ELEMENTS AYANT UNE FACE COMMUNE
! |                |    | AVEC L'ELEMENT .  SI IFABOR<0 OU NUL         |
! |                |    | ON A UNE FACE LIQUIDE,SOLIDE,OU PERIODIQUE   |
! |  X..,Y..,ZPLOT |<-->| POSITIONS SUCCESSIVES DES DERIVANTS.         |
! |    DX,DY,DZ    | -- | STOCKAGE DES SOUS-PAS .                      |
! |    SHP         |<-->| COORDONNEES BARYCENTRIQUES 2D AU PIED DES    |
! |                |    | COURBES CARACTERISTIQUES.                    |
! |    SHZ         |<-->| COORDONNEES BARYCENTRIQUES SUIVANT Z DES     |
! |                |    | NOEUDS DANS LEURS ETAGES "ETA" ASSOCIES.     |
! |    ELT         |<-->| NUMEROS DES ELEMENTS 2D CHOISIS POUR CHAQUE  |
! |                |    | NOEUD.                                       |
! |    ETA         |<-->| NUMEROS DES ETAGES CHOISIS POUR CHAQUE NOEUD.|
! |    NPLOT       | -->| NOMBRE DE DERIVANTS.                         |
! |    NPOIN2      | -->| NOMBRE DE POINTS DU MAILLAGE 2D.             |
! |    NELEM2      | -->| NUMBER OF ELEMENTS IN THE 2D MESH
! |    NELMAX2     | -->| MAXIMUM NUMBER OF ELEMENTS IN THE 2D MESH
! |    NPLAN       | -->| NOMBRE DE PLANS.                             |
! |    SURDET      | -->| VARIABLE UTILISEE PAR LA TRANSFORMEE ISOPARAM.
! |    SENS        | -->| DESCENTE OU REMONTEE DES CARACTERISTIQUES.   |
! |    ISO         | -- | STOCKAGE BINAIRE DE LA FACE DE SORTIE.       |
! |________________|____|______________________________________________|
!  MODE: -->(DONNEE NON MODIFIEE),<--(RESULTAT),<-->(DONNEE MODIFIEE)
!-----------------------------------------------------------------------
!     - APPELE PAR : CARACT
!     - PROGRAMMES APPELES : NEANT
!
!***********************************************************************
!
      USE bief
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)    :: SENS,NPLAN,NCHDIM,NELMAX2
      INTEGER         , INTENT(IN)    :: NPOIN2,NELEM2,NPLOT,NRK
      INTEGER         , INTENT(IN)    :: IKLE2(nelmax2,3)
      INTEGER         , INTENT(INOUT) :: ELT(nplot),NCHARA
      DOUBLE PRECISION, INTENT(IN)    :: U(npoin2,nplan),V(npoin2,nplan)
      DOUBLE PRECISION, INTENT(IN)    :: W(npoin2,nplan),SURDET(nelem2)
      DOUBLE PRECISION, INTENT(INOUT) :: XPLOT(nplot),YPLOT(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: ZPLOT(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: SHP(3,nplot),SHZ(nplot)
      DOUBLE PRECISION, INTENT(IN)    :: X(npoin2),Y(npoin2),DT
      DOUBLE PRECISION, INTENT(IN)    :: Z(npoin2,nplan),ZSTAR(nplan)
      DOUBLE PRECISION, INTENT(INOUT) :: DX(nplot),DY(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: DZ(nplot)
      INTEGER         , INTENT(IN)    :: IBOR(nelmax2,5,nplan-1)
      INTEGER         , INTENT(INOUT) :: ETA(nplot)
      INTEGER         , INTENT(IN)    :: IFAPAR(6,*)
      LOGICAL, INTENT(IN)             :: ADD
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELE,ISO,ISPDONE,NSP
      INTEGER :: IPLOT,ISP,I1,I2,I3,IEL,IET,IET2,ISOH,ISOV,IFA,ISUI(3)
!
      DOUBLE PRECISION PAS,A1,DX1,DY1,DXP,DYP,XP,YP,ZP,DENOM
      DOUBLE PRECISION DELTAZ,PAS2
!
      INTRINSIC abs , int , max , sqrt
!
      parameter( isui = (/ 2 , 3 , 1 /) )
      DOUBLE PRECISION, PARAMETER :: EPSILO = -1.d-6
      DOUBLE PRECISION, PARAMETER :: EPSDZ = 1.d-4
!
!-----------------------------------------------------------------------
!     FOR EVERY POINT
!-----------------------------------------------------------------------
!
      DO iplot = 1 , nplot
!
        IF(add) THEN
!
        xplot(iplot)   = recvchar(iplot)%XP
        yplot(iplot)   = recvchar(iplot)%YP
        zplot(iplot)   = recvchar(iplot)%ZP
        dx(iplot)      = recvchar(iplot)%DX
        dy(iplot)      = recvchar(iplot)%DY
        dz(iplot)      = recvchar(iplot)%DZ
        elt(iplot)     = recvchar(iplot)%INE
        eta(iplot)     = recvchar(iplot)%KNE
        nsp            = recvchar(iplot)%NSP ! R-K STEPS TO BE FULLFILLED
        ispdone        = recvchar(iplot)%ISP ! R-K STEPS ALREADY DONE
        iel = elt(iplot)
        iet = eta(iplot)
        xp  = xplot(iplot)
        yp  = yplot(iplot)
        zp  = zplot(iplot)
        i1 = ikle2(iel,1)
        i2 = ikle2(iel,2)
        i3 = ikle2(iel,3)
        shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                 -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
        shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                 -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
        shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                 -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
        shz(iplot) = (zp-zstar(iet)) / (zstar(iet+1)-zstar(iet))
!       ASSUME ALL ARE LOCALISED, IT WILL BE SET OTHERWISE IF LOST-AGAIN
        recvchar(iplot)%NEPID=-1
!
        ELSE
!
          iel = elt(iplot)
!         POINTS WITH IEL=0 ARE TREATED SO THAT THE FINAL
!         INTERPOLATION GIVES 0.,
!         AND WE SKIP TO NEXT POINT IPLOT (CYCLE)
          IF(iel.EQ.0) THEN
            elt(iplot)=1
            eta(iplot)=1
            shp(1,iplot)=0.d0
            shp(2,iplot)=0.d0
            shp(3,iplot)=0.d0
            shz(iplot)=0.d0
            cycle
          ENDIF
          iet = eta(iplot)
          i1 = ikle2(iel,1)
          i2 = ikle2(iel,2)
          i3 = ikle2(iel,3)
          dxp =( u(i1,iet  )*shp(1,iplot)
     &         + u(i2,iet  )*shp(2,iplot)
     &         + u(i3,iet  )*shp(3,iplot) )*(1.d0-shz(iplot))
     &        +( u(i1,iet+1)*shp(1,iplot)
     &         + u(i2,iet+1)*shp(2,iplot)
     &         + u(i3,iet+1)*shp(3,iplot) )*shz(iplot)
          dyp =( v(i1,iet  )*shp(1,iplot)
     &         + v(i2,iet  )*shp(2,iplot)
     &         + v(i3,iet  )*shp(3,iplot) )*(1.d0-shz(iplot))
     &        +( v(i1,iet+1)*shp(1,iplot)
     &         + v(i2,iet+1)*shp(2,iplot)
     &         + v(i3,iet+1)*shp(3,iplot) )*shz(iplot)
          nsp=max(1,int(nrk*dt*sqrt((dxp**2+dyp**2)*surdet(iel))))
          ispdone=1
!
        ENDIF
!
        pas = sens * dt / nsp
!
!       LOOP ON RUNGE-KUTTA SUB-STEPS
!
!       COMPILER MUST DO NOTHING IF ISPDONE>NSP
!       IN MODE "ADD", ISP = ISPDONE HAS NOT BEEN FULLY DONE
!       IT IS RESTARTED HERE
!
      DO isp = ispdone,nsp
!
!-----------------------------------------------------------------------
!       LOCALISING THE ARRIVAL POINT
!-----------------------------------------------------------------------
!
        pas2 = pas
!
!       IN MODE "ADD" ITERATIONS ALREADY DONE ARE SKIPPED AND
!                     CHARACTERISTICS GONE IN ANOTHER SUB-DOMAIN SKIPPED
!
        IF(add) THEN
          IF(isp.EQ.ispdone) GO TO 50
          IF(recvchar(iplot)%NEPID.NE.-1) cycle
        ENDIF
!
        iel = elt(iplot)
        iet = eta(iplot)
        i1 = ikle2(iel,1)
        i2 = ikle2(iel,2)
        i3 = ikle2(iel,3)
!
        dx(iplot) = ((u(i1,iet  )*shp(1,iplot)
     &              + u(i2,iet  )*shp(2,iplot)
     &              + u(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &              +(u(i1,iet+1)*shp(1,iplot)
     &              + u(i2,iet+1)*shp(2,iplot)
     &              + u(i3,iet+1)*shp(3,iplot))*shz(iplot) ) * pas
!
        dy(iplot) = ((v(i1,iet  )*shp(1,iplot)
     &              + v(i2,iet  )*shp(2,iplot)
     &              + v(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &              +(v(i1,iet+1)*shp(1,iplot)
     &              + v(i2,iet+1)*shp(2,iplot)
     &              + v(i3,iet+1)*shp(3,iplot))*shz(iplot) ) * pas
!
        deltaz =  (z(i1,iet+1)-z(i1,iet))*shp(1,iplot)
     &          + (z(i2,iet+1)-z(i2,iet))*shp(2,iplot)
     &          + (z(i3,iet+1)-z(i3,iet))*shp(3,iplot)
!
        IF(deltaz.GT.epsdz) THEN
!         DIVISION BY DELTAZ IS DUE TO THE FACT THAT W IS
!         W* MULTIPLIED BY DELTAZ (IT STEMS FROM TRIDW2 IN TELEMAC3D)
          dz(iplot) = ((w(i1,iet  )*shp(1,iplot)
     &                + w(i2,iet  )*shp(2,iplot)
     &                + w(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &                +(w(i1,iet+1)*shp(1,iplot)
     &                + w(i2,iet+1)*shp(2,iplot)
     &                + w(i3,iet+1)*shp(3,iplot))*shz(iplot) )
     &                * pas * (zstar(iet+1)-zstar(iet)) / deltaz
        ELSE
          dz(iplot) = 0.d0
        ENDIF
!
        xp = xplot(iplot) + dx(iplot)
        yp = yplot(iplot) + dy(iplot)
        zp = zplot(iplot) + dz(iplot)
!
        shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                 -(y(i3)-y(i2))*(xp-x(i2))) * surdet(iel)
        shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                 -(y(i1)-y(i3))*(xp-x(i3))) * surdet(iel)
        shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                 -(y(i2)-y(i1))*(xp-x(i1))) * surdet(iel)
!
        shz(iplot) = (zp-zstar(iet)) / (zstar(iet+1)-zstar(iet))
!
        xplot(iplot) = xp
        yplot(iplot) = yp
        zplot(iplot) = zp
!
        IF(add) THEN
!         CONTINUOUS SETTING OF THE REACHED POSITION FOR IPLOT
!         AND THE NUMBER OF STEPS DONE ALREADY
          recvchar(iplot)%XP=xplot(iplot)
          recvchar(iplot)%YP=yplot(iplot)
          recvchar(iplot)%ZP=zplot(iplot)
          recvchar(iplot)%DX=dx(iplot)
          recvchar(iplot)%DY=dy(iplot)
          recvchar(iplot)%DZ=dz(iplot)
          recvchar(iplot)%INE=elt(iplot)
          recvchar(iplot)%ISP=isp
        ENDIF
!
!-----------------------------------------------------------------------
!       TEST: IS THE PATHLINE WENT OUT THE ORIGINAL ELEMENT
!-----------------------------------------------------------------------
!
50      CONTINUE
!
        iso = 0
        IF(shp(1,iplot).LT.epsilo)      iso=ibset(iso,2)
        IF(shp(2,iplot).LT.epsilo)      iso=ibset(iso,3)
        IF(shp(3,iplot).LT.epsilo)      iso=ibset(iso,4)
        IF(shz(iplot)  .LT.epsilo)      iso=ibset(iso,0)
        IF(shz(iplot)  .GT.1.d0-epsilo) iso=ibset(iso,1)
!
        IF(iso.NE.0) THEN
!
!-----------------------------------------------------------------------
!         HERE WE ARE OUT OF THE ELEMENT
!-----------------------------------------------------------------------
!
          isoh = iand(iso,28)
          isov = iand(iso, 3)
          iel = elt(iplot)
          iet = eta(iplot)
          xp = xplot(iplot)
          yp = yplot(iplot)
          zp = zplot(iplot)
!
          IF(isoh.NE.0) THEN
!
            IF(isoh.EQ.4) THEN
              ifa = 2
            ELSEIF(isoh.EQ.8) THEN
              ifa = 3
            ELSEIF(isoh.EQ.16) THEN
              ifa = 1
            ELSEIF(isoh.EQ.12) THEN
              ifa = 2
              IF(dx(iplot)*(y(ikle2(iel,3))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,3))-xp)) ifa = 3
            ELSEIF(isoh.EQ.24) THEN
              ifa = 3
              IF(dx(iplot)*(y(ikle2(iel,1))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,1))-xp)) ifa = 1
            ELSE
              ifa = 1
              IF(dx(iplot)*(y(ikle2(iel,2))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,2))-xp)) ifa = 2
            ENDIF
!
            IF(isov.GT.0) THEN
              IF(abs(dz(iplot)).GT.epsdz) THEN
                a1 = (zp-zstar(iet+isov-1)) / dz(iplot)
              ELSE
                a1 = 0.d0
              ENDIF
              i1 = ikle2(iel,ifa)
              i2 = ikle2(iel,isui(ifa))
              IF ((x(i2)-x(i1))*(yp-a1*dy(iplot)-y(i1)).GT.
     &            (y(i2)-y(i1))*(xp-a1*dx(iplot)-x(i1))) ifa=isov+3
            ENDIF
!
          ELSE
!
            ifa = isov + 3
!
          ENDIF
!
          iel = ibor(iel,ifa,iet)
!
          IF(ifa.LE.3) THEN
!
!-----------------------------------------------------------------------
!           HERE WE ARRIVE IN ANOTHER ELEMENT THROUGH A QUADRANGULAR FACE
!-----------------------------------------------------------------------
!
            IF(iel.GT.0) THEN
!
!-----------------------------------------------------------------------
!             RELOCALISING IN ADJACENT ELEMENT
!-----------------------------------------------------------------------
!
              i1 = ikle2(iel,1)
              i2 = ikle2(iel,2)
              i3 = ikle2(iel,3)
              elt(iplot) = iel
              shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
              shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
              shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
!
              GOTO 50
!
            ENDIF
!
!-----------------------------------------------------------------------
!           HERE WE PASS TO A NEIGHBOUR SUBDOMAIN AND COLLECT DATA
!-----------------------------------------------------------------------
!
          IF(iel.EQ.-2) THEN
            IF(.NOT.add) THEN
!             INTERFACE CROSSING
              CALL collect_char
     &            (ipid,iplot,elt(iplot),ifa,eta(iplot),0,isp,
     &             nsp,xplot(iplot),yplot(iplot),
     &             zplot(iplot),0.d0,
     &             dx(iplot),dy(iplot),dz(iplot),0.d0,
     &             ifapar,nchdim,nchara)
            ELSE
!             A LOST-AGAIN TRACEBACK DETECTED
!             PROCESSOR NUMBER
              recvchar(iplot)%NEPID=ifapar(ifa,elt(iplot))
              recvchar(iplot)%INE=ifapar(ifa+3,elt(iplot))
              recvchar(iplot)%KNE=eta(iplot)
            ENDIF
!           EXITING LOOP ON ISP
            EXIT
          ENDIF
!
!-----------------------------------------------------------------------
!           SPECIAL TREATMENT FOR SOLID OR LIQUID BOUNDARIES
!-----------------------------------------------------------------------
!
            dxp = dx(iplot)
            dyp = dy(iplot)
            i1  = ikle2(elt(iplot),ifa)
            i2  = ikle2(elt(iplot),isui(ifa))
            dx1 = x(i2) - x(i1)
            dy1 = y(i2) - y(i1)
!
            IF(iel.EQ.-1) THEN
!
!-----------------------------------------------------------------------
!             HERE SOLID BOUNDARY, VELOCITY IS PROJECTED ON THE BOUNDARY
!             AND WE GO ON
!-----------------------------------------------------------------------
!
              a1 = (dxp*dx1+dyp*dy1) / (dx1**2+dy1**2)
              dx(iplot) = a1 * dx1
              dy(iplot) = a1 * dy1
!
              a1=((xp-x(i1))*dx1+(yp-y(i1))*dy1)/(dx1**2+dy1**2)
              shp(          ifa  ,iplot) = 1.d0 - a1
              shp(     isui(ifa) ,iplot) = a1
              shp(isui(isui(ifa)),iplot) = 0.d0
              xplot(iplot) = x(i1) + a1 * dx1
              yplot(iplot) = y(i1) + a1 * dy1
              IF(add) THEN
                recvchar(iplot)%XP=xplot(iplot)
                recvchar(iplot)%YP=yplot(iplot)
                recvchar(iplot)%ZP=zplot(iplot)
                recvchar(iplot)%DX=dx(iplot)
                recvchar(iplot)%DY=dy(iplot)
                recvchar(iplot)%DZ=dz(iplot)
              ENDIF
!
              GOTO 50
!
            ELSEIF(iel.EQ.0) THEN
!
!-----------------------------------------------------------------------
!           HERE WE HAVE A LIQUID BOUNDARY, THE CHARACTERISTIC IS STOPPED
!-----------------------------------------------------------------------
!
            denom = dxp*dy1-dyp*dx1
            IF(abs(denom).GT.1.d-12) THEN
              a1 = (dxp*(yp-y(i1))-dyp*(xp-x(i1))) / denom
            ELSE
              a1 = 0.d0
            ENDIF
            IF(a1.GT.1.d0) a1 = 1.d0
            IF(a1.LT.0.d0) a1 = 0.d0
            shp(          ifa  ,iplot) = 1.d0 - a1
            shp(     isui(ifa) ,iplot) = a1
            shp(isui(isui(ifa)),iplot) = 0.d0
            xplot(iplot) = x(i1) + a1 * dx1
            yplot(iplot) = y(i1) + a1 * dy1
            IF(abs(dxp).GT.abs(dyp)) THEN
              a1 = (xp-xplot(iplot))/dxp
            ELSE
              a1 = (yp-yplot(iplot))/dyp
            ENDIF
            zplot(iplot) = zp - a1*dz(iplot)
            shz(iplot) = (zplot(iplot)-zstar(iet))
     &                 / (zstar(iet+1)-zstar(iet))
!           THIS IS A MARKER FOR PARTICLES EXITING A DOMAIN
!           SENS=-1 FOR BACKWARD CHARACTERISTICS
            elt(iplot) = - sens * elt(iplot)
!           EXITING LOOP ON ISP
            EXIT
!
            ELSE
!
              WRITE(lu,*) 'UNEXPECTED CASE IN SCHAR41_SIGMA'
              WRITE(lu,*) 'IEL=',iel
              CALL plante(1)
              stop
!
            ENDIF
!
          ELSE
!
!-----------------------------------------------------------------------
!  CASE IFA = 4 OR 5
!  HERE WE EXIT THROUGH TOP OR BOTTOM OF THE PRISM
!-----------------------------------------------------------------------
!
            ifa = ifa - 4
!           HENCE IFA NOW EQUALS 0 OR 1
!
            IF(iel.EQ.1) THEN
!
!-----------------------------------------------------------------------
!           NO NEED TO RECOMPUTE THE VELOCITIES,
!           RELOCALISING IN NEW ELEMENT
!-----------------------------------------------------------------------
!
              eta(iplot) = iet + ifa + ifa - 1
              shz(iplot) = (zp-zstar(eta(iplot)))
     &                   / (zstar(eta(iplot)+1)-zstar(eta(iplot)))
!
              IF(add) THEN
                recvchar(iplot)%KNE=eta(iplot)
              ENDIF
!
              iso = isoh
!
              IF(shz(iplot).LT.     epsilo) iso=ibset(iso,0)
              IF(shz(iplot).GT.1.d0-epsilo) iso=ibset(iso,1)
!
              GOTO 50
!
            ENDIF
!
            IF(iel.EQ.-1) THEN
!
!-----------------------------------------------------------------------
!  LA, ON SAIT QUE LA FACE DE SORTIE EST UNE FRONTIERE SOLIDE
!  ON PROJETTE LE RELIQUAT SUR LA FRONTIERE PUIS ON SE RELOCALISE
!-----------------------------------------------------------------------
!
              zplot(iplot) = zstar(iet+ifa)
              dz(iplot) = 0.d0
              shz(iplot) = ifa
!
              iso = isoh
              IF(isoh.NE.0) GOTO 50
!
            ELSE
!
!-----------------------------------------------------------------------
!  LA, ON SAIT QUE LA FACE DE SORTIE EST UNE FRONTIERE LIQUIDE (CAS 0)
!      ON ARRETE ALORS LA REMONTEE DES CARACTERISTIQUES (SIGNE DE ELT)
!  OU, QUE L'ON VIENT DE TRAVERSER UN PLAN AVEC RECALCUL DES VITESSES
!  DEMANDE (CAS 2)
!-----------------------------------------------------------------------
!
              IF(abs(dz(iplot)).GT.epsdz) THEN
                a1 = (zp-zstar(iet+ifa)) / dz(iplot)
              ELSE
                a1 = 0.d0
              ENDIF
              xp = xp - a1*dx(iplot)
              yp = yp - a1*dy(iplot)
              zp = zstar(iet+ifa)
              iele = elt(iplot)
              i1 = ikle2(iele,1)
              i2 = ikle2(iele,2)
              i3 = ikle2(iele,3)
!
              shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iele)
              shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iele)
              shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iele)
!
              IF(iel.EQ.2) THEN
!
!-----------------------------------------------------------------------
!  LA, ON SAIT QUE LA FACE DE SORTIE SE SITUE SUR LE PLAN OU ON DEMANDE
!  UN RECALCUL DES VITESSES
!-----------------------------------------------------------------------
!
!               IF IFA = 1 EXIT THROUGH THE TOP
!               IF IFA = 0 EXIT THROUGH THE BOTTOM
!               THEN NEW IET IS  IET+1 IF IFA = 1
!                            AND IET-1 IF IFA = 0
!               THIS IS SUMMARISED BY IET=IET+2*IFA-1
!
!               RECOMPUTED VELOCITIES MUST BE TAKEN AT IET2=IET+IFA
!               I.E. BOTTOM IF EXIT THROUGH THE BOTTOM
!               AND TOP IF EXIT THROUGH THE TOP
!
                iet2 = iet + ifa
                iet  = iet + ifa + ifa - 1
!
                pas2 = pas2 * a1
!
                dx(iplot) = ( u(i1,iet2)*shp(1,iplot)
     &                      + u(i2,iet2)*shp(2,iplot)
     &                      + u(i3,iet2)*shp(3,iplot) ) * pas2
!
                dy(iplot) = ( v(i1,iet2)*shp(1,iplot)
     &                      + v(i2,iet2)*shp(2,iplot)
     &                      + v(i3,iet2)*shp(3,iplot) ) * pas2
!
                deltaz =  (z(i1,iet+1)-z(i1,iet))*shp(1,iplot)
     &                  + (z(i2,iet+1)-z(i2,iet))*shp(2,iplot)
     &                  + (z(i3,iet+1)-z(i3,iet))*shp(3,iplot)
!
                IF(deltaz.GT.epsdz) THEN
                  dz(iplot) = ( w(i1,iet2)*shp(1,iplot)
     &                        + w(i2,iet2)*shp(2,iplot)
     &                        + w(i3,iet2)*shp(3,iplot) ) * pas2
     &                        * (zstar(iet+1)-zstar(iet)) / deltaz
                ELSE
                  dz(iplot) = 0.d0
                ENDIF
!
                xp = xp + dx(iplot)
                yp = yp + dy(iplot)
                zp = zp + dz(iplot)
!
                shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2))) * surdet(iele)
                shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3))) * surdet(iele)
                shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1))) * surdet(iele)
                shz(iplot)=(zp-zstar(iet))/(zstar(iet+1)-zstar(iet))
!
                xplot(iplot) = xp
                yplot(iplot) = yp
                zplot(iplot) = zp
                eta(iplot) = iet
!
                iso = 0
!
                IF(shp(1,iplot).LT.epsilo) iso=ibset(iso,2)
                IF(shp(2,iplot).LT.epsilo) iso=ibset(iso,3)
                IF(shp(3,iplot).LT.epsilo) iso=ibset(iso,4)
!
                IF(shz(iplot).LT.     epsilo) iso=ibset(iso,0)
                IF(shz(iplot).GT.1.d0-epsilo) iso=ibset(iso,1)
!
                IF(add) THEN
                  recvchar(iplot)%XP=xplot(iplot)
                  recvchar(iplot)%YP=yplot(iplot)
                  recvchar(iplot)%ZP=zplot(iplot)
                  recvchar(iplot)%DX=dx(iplot)
                  recvchar(iplot)%DY=dy(iplot)
                  recvchar(iplot)%DZ=dz(iplot)
                ENDIF
!
                GOTO 50
!
              ENDIF
!
              xplot(iplot) = xp
              yplot(iplot) = yp
              zplot(iplot) = zp
              shz(iplot) = ifa
!             THIS IS A MARKER FOR PARTICLES EXITING A DOMAIN
!             SENS=-1 FOR BACKWARD CHARACTERISTICS
              elt(iplot) = - sens * elt(iplot)
!             EXITING LOOP ON ISP
              EXIT
!
            ENDIF
!
          ENDIF
!
!       IF(ISO.NE.0)
        ENDIF
!
        ENDDO
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE schar41_sigma
!                       **********************
                        SUBROUTINE schar41_per
!                       **********************
!
     &( u , v , w , dt , nrk , x , y , zstar , ikle2 , ibor ,
     &  xplot , yplot , zplot , dx , dy , dz , shp , shz , elt , eta ,
     &  nplot , npoin2 , nelem2 , nelmax2,nplan , surdet ,
     &  sens  , ifapar, nchdim,nchara,add)
!
!***********************************************************************
! BIEF VERSION 7.3
!
!brief    LIKE SCHAR41 BUT WITH PERIODICITY ON THE VERTICAL
!         DIFFERENCES MAEKED WITH "PERIODICITY"
!
!
!history  J-M HERVOUET (LNHE)
!+        31/07/2012
!+        V6P3
!+      First version (differences taken in Tomawac)
!
!***********************************************************************
!
!  FONCTION :
!
!     REMONTEE OU DESCENTE
!     DES COURBES CARACTERISTIQUES
!     SUR DES PRISMES DE TELEMAC-3D
!     DANS L'INTERVALLE DE TEMPS DT
!     AVEC UNE DISCRETISATION ELEMENTS FINIS
!
!
!  DISCRETISATION :
!
!     LE DOMAINE EST APPROCHE PAR UNE DISCRETISATION ELEMENTS FINIS
!     UNE APPROXIMATION LOCALE EST DEFINIE POUR LE VECTEUR VITESSE :
!     LA VALEUR EN UN POINT D'UN ELEMENT NE DEPEND QUE DES VALEURS
!     AUX NOEUDS DE CET ELEMENT
!
!
!  RESTRICTIONS ET HYPOTHESES :
!
!     LE CHAMP CONVECTEUR U EST SUPPOSE INDEPENDANT DU TEMPS
!
!-----------------------------------------------------------------------
!                             ARGUMENTS
! .________________.____.______________________________________________.
! |      NOM       |MODE|                   ROLE                       |
! |________________|____|______________________________________________|
! |    U,V,W       | -->| COMPONENTS OF ADVECTION VELOCITY
! |                |    | BUT W IS W* x DELTAZ (STEMS FROM TRIDW2)
! |    DT          | -->| PAS DE TEMPS.                                |
! |    NRK         | -->| NOMBRE DE SOUS-PAS DE RUNGE-KUTTA.           |
! |    X,Y,ZSTAR   | -->| COORDONNEES DES POINTS DU MAILLAGE.
! |    IKLE2       | -->| TRANSITION ENTRE LES NUMEROTATIONS LOCALE    |
! |                |    | ET GLOBALE DU MAILLAGE 2D.                   |
! |    IBOR        | -->| NUMEROS 2D DES ELEMENTS AYANT UNE FACE COMMUNE
! |                |    | AVEC L'ELEMENT .  SI IFABOR<0 OU NUL         |
! |                |    | ON A UNE FACE LIQUIDE,SOLIDE,OU PERIODIQUE   |
! |  X..,Y..,ZPLOT |<-->| POSITIONS SUCCESSIVES DES DERIVANTS.         |
! |    DX,DY,DZ    | -- | STOCKAGE DES SOUS-PAS .                      |
! |    SHP         |<-->| COORDONNEES BARYCENTRIQUES 2D AU PIED DES    |
! |                |    | COURBES CARACTERISTIQUES.                    |
! |    SHZ         |<-->| COORDONNEES BARYCENTRIQUES SUIVANT Z DES     |
! |                |    | NOEUDS DANS LEURS ETAGES "ETA" ASSOCIES.     |
! |    ELT         |<-->| NUMEROS DES ELEMENTS 2D CHOISIS POUR CHAQUE  |
! |                |    | NOEUD.                                       |
! |    ETA         |<-->| NUMEROS DES ETAGES CHOISIS POUR CHAQUE NOEUD.|
! |    NPLOT       | -->| NOMBRE DE DERIVANTS.                         |
! |    NPOIN2      | -->| NOMBRE DE POINTS DU MAILLAGE 2D.             |
! |    NELEM2      | -->| NOMBRE D'ELEMENTS DU MAILLAGE 2D.            |
! |    NPLAN       | -->| NOMBRE DE PLANS.                             |
! |    SURDET      | -->| VARIABLE UTILISEE PAR LA TRANSFORMEE ISOPARAM.
! |    SENS        | -->| DESCENTE OU REMONTEE DES CARACTERISTIQUES.   |
! |    ISO         | -- | STOCKAGE BINAIRE DE LA FACE DE SORTIE.       |
! |________________|____|______________________________________________|
!  MODE: -->(DONNEE NON MODIFIEE),<--(RESULTAT),<-->(DONNEE MODIFIEE)
!-----------------------------------------------------------------------
!     - APPELE PAR : CARACT
!     - PROGRAMMES APPELES : NEANT
!
!***********************************************************************
!
      USE bief
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)    :: SENS,NPLAN,NCHDIM,NELMAX2
      INTEGER         , INTENT(IN)    :: NPOIN2,NELEM2,NPLOT,NRK
      INTEGER         , INTENT(IN)    :: IKLE2(nelmax2,3)
      INTEGER         , INTENT(INOUT) :: ELT(nplot),NCHARA
      DOUBLE PRECISION, INTENT(IN)    :: U(npoin2,nplan),V(npoin2,nplan)
      DOUBLE PRECISION, INTENT(IN)    :: W(npoin2,nplan),SURDET(nelem2)
      DOUBLE PRECISION, INTENT(INOUT) :: XPLOT(nplot),YPLOT(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: ZPLOT(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: SHP(3,nplot),SHZ(nplot)
      DOUBLE PRECISION, INTENT(IN)    :: X(npoin2),Y(npoin2),DT
!                                              PERIODICITY
      DOUBLE PRECISION, INTENT(IN)    :: ZSTAR(nplan+1)
      DOUBLE PRECISION, INTENT(INOUT) :: DX(nplot),DY(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: DZ(nplot)
      INTEGER         , INTENT(IN)    :: IBOR(nelmax2,5,nplan-1)
      INTEGER         , INTENT(INOUT) :: ETA(nplot)
      INTEGER         , INTENT(IN)    :: IFAPAR(6,*)
      LOGICAL         , INTENT(IN)    :: ADD
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELE,ISO,ISPDONE,NSP,NSPMAX,IETP1
      INTEGER IPLOT,ISP,I1,I2,I3,IEL,IET,IET2,ISOH,ISOV,IFA,ISUI(3)
!
      DOUBLE PRECISION PAS,A1,DX1,DY1,DXP,DYP,DZP,XP,YP,ZP,DENOM,PAS2
!
      INTRINSIC abs , int , max , sqrt
!
      parameter( isui = (/ 2 , 3 , 1 /) )
      DOUBLE PRECISION, PARAMETER :: EPSILO = -1.d-6
      DOUBLE PRECISION, PARAMETER :: EPSDZ = 1.d-4
!
      nspmax=1
      eta1(nplan)=1
!
!-----------------------------------------------------------------------
!     FOR EVERY POINT
!-----------------------------------------------------------------------
!
      DO iplot = 1 , nplot
!
        IF(add) THEN
!
        xplot(iplot)   = recvchar(iplot)%XP
        yplot(iplot)   = recvchar(iplot)%YP
        zplot(iplot)   = recvchar(iplot)%ZP
        dx(iplot)      = recvchar(iplot)%DX
        dy(iplot)      = recvchar(iplot)%DY
        dz(iplot)      = recvchar(iplot)%DZ
        elt(iplot)     = recvchar(iplot)%INE
        eta(iplot)     = recvchar(iplot)%KNE
        nsp            = recvchar(iplot)%NSP ! R-K STEPS TO BE FULLFILLED
        ispdone        = recvchar(iplot)%ISP ! R-K STEPS ALREADY DONE
        iel = elt(iplot)
        iet = eta(iplot)
        xp  = xplot(iplot)
        yp  = yplot(iplot)
        zp  = zplot(iplot)
        i1 = ikle2(iel,1)
        i2 = ikle2(iel,2)
        i3 = ikle2(iel,3)
        shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                 -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
        shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                 -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
        shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                 -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
        shz(iplot) = (zp-zstar(iet)) / (zstar(iet+1)-zstar(iet))
!       ASSUME ALL ARE LOCALISED, IT WILL BE SET OTHERWISE IF LOST-AGAIN
        recvchar(iplot)%NEPID=-1
!
        ELSE
!
          iel = elt(iplot)
!         POINTS WITH IEL=0 ARE TREATED SO THAT THE FINAL
!         INTERPOLATION GIVES 0.,
!         AND WE SKIP TO NEXT POINT IPLOT (CYCLE)
          IF(iel.EQ.0) THEN
            elt(iplot)=1
            eta(iplot)=1
            shp(1,iplot)=0.d0
            shp(2,iplot)=0.d0
            shp(3,iplot)=0.d0
            shz(iplot)=0.d0
            cycle
          ENDIF
          iet = eta(iplot)
          i1 = ikle2(iel,1)
          i2 = ikle2(iel,2)
          i3 = ikle2(iel,3)
!         HERE IET+1 IS ALWAYS < NPLAN+1 (SEE GTSH41)
          dxp = u(i1,iet  )*shp(1,iplot)*(1.d0-shz(iplot))
     &        + u(i2,iet  )*shp(2,iplot)*(1.d0-shz(iplot))
     &        + u(i3,iet  )*shp(3,iplot)*(1.d0-shz(iplot))
     &        + u(i1,iet+1)*shp(1,iplot)*shz(iplot)
     &        + u(i2,iet+1)*shp(2,iplot)*shz(iplot)
     &        + u(i3,iet+1)*shp(3,iplot)*shz(iplot)
          dyp = v(i1,iet  )*shp(1,iplot)*(1.d0-shz(iplot))
     &        + v(i2,iet  )*shp(2,iplot)*(1.d0-shz(iplot))
     &        + v(i3,iet  )*shp(3,iplot)*(1.d0-shz(iplot))
     &        + v(i1,iet+1)*shp(1,iplot)*shz(iplot)
     &        + v(i2,iet+1)*shp(2,iplot)*shz(iplot)
     &        + v(i3,iet+1)*shp(3,iplot)*shz(iplot)
          nsp=max(1,int(nrk*dt*sqrt((dxp**2+dyp**2)*surdet(iel))))
!!        PERIODICITY
          dzp = w(i1,iet  )*shp(1,iplot)*(1.d0-shz(iplot))
     &        + w(i2,iet  )*shp(2,iplot)*(1.d0-shz(iplot))
     &        + w(i3,iet  )*shp(3,iplot)*(1.d0-shz(iplot))
     &        + w(i1,iet+1)*shp(1,iplot)*shz(iplot)
     &        + w(i2,iet+1)*shp(2,iplot)*shz(iplot)
     &        + w(i3,iet+1)*shp(3,iplot)*shz(iplot)
!
          nsp=max(nsp,int(nrk*dt*abs(dzp/(zstar(iet+1)-zstar(iet)))))
          nspmax = max(nspmax,nsp)
!!        END PERIODICITY
!
          ispdone=1
!
        ENDIF
!
        pas = sens * dt / nsp
!
!       LOOP ON RUNGE-KUTTA SUB-STEPS
!
!       COMPILER MUST DO NOTHING IF ISPDONE>NSP
!       IN MODE "ADD", ISP = ISPDONE HAS NOT BEEN FULLY DONE
!       IT IS RESTARTED HERE
!
      DO isp = ispdone,nsp
!
!-----------------------------------------------------------------------
!       LOCALISING THE ARRIVAL POINT
!-----------------------------------------------------------------------
!
        pas2 = pas
!
!       IN MODE "ADD" ITERATIONS ALREADY DONE ARE SKIPPED AND
!                     CHARACTERISTICS GONE IN ANOTHER SUB-DOMAIN SKIPPED
!
        IF(add) THEN
          IF(isp.EQ.ispdone) GO TO 50
          IF(recvchar(iplot)%NEPID.NE.-1) cycle
        ENDIF
!
        iel = elt(iplot)
        iet = eta(iplot)
!!      PERIODICITY IETP1 REPLACES IET+1 (BUT NOT ALWAYS)
        ietp1=eta1(iet)
!
        i1 = ikle2(iel,1)
        i2 = ikle2(iel,2)
        i3 = ikle2(iel,3)
!
        dx(iplot) = ((u(i1,iet  )*shp(1,iplot)
     &              + u(i2,iet  )*shp(2,iplot)
     &              + u(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &              +(u(i1,ietp1)*shp(1,iplot)
     &              + u(i2,ietp1)*shp(2,iplot)
     &              + u(i3,ietp1)*shp(3,iplot))*shz(iplot) )*pas
!
        dy(iplot) = ((v(i1,iet  )*shp(1,iplot)
     &              + v(i2,iet  )*shp(2,iplot)
     &              + v(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &              +(v(i1,ietp1)*shp(1,iplot)
     &              + v(i2,ietp1)*shp(2,iplot)
     &              + v(i3,ietp1)*shp(3,iplot))*shz(iplot) )*pas
!! PERIODICITY, NO DELTA, Z NEVER USED
!       DELTAZ =  (Z(I1,IET+1)-Z(I1,IET))*SHP(1,IPLOT)
!    &          + (Z(I2,IET+1)-Z(I2,IET))*SHP(2,IPLOT)
!    &          + (Z(I3,IET+1)-Z(I3,IET))*SHP(3,IPLOT)
!
!       IF(DELTAZ.GT.EPSDZ) THEN
!         DIVISION BY DELTAZ IS DUE TO THE FACT THAT W IS
!         W* MULTIPLIED BY DELTAZ (IT STEMS FROM TRIDW2 IN TELEMAC3D)
          dz(iplot) = ((w(i1,iet  )*shp(1,iplot)
     &                + w(i2,iet  )*shp(2,iplot)
     &                + w(i3,iet  )*shp(3,iplot))*(1.d0-shz(iplot))
     &                +(w(i1,ietp1)*shp(1,iplot)
     &                + w(i2,ietp1)*shp(2,iplot)
     &                + w(i3,ietp1)*shp(3,iplot))*shz(iplot))
!! PERIODICITY, PAS DE DELTAZ
!    &                * PAS * (ZSTAR(IET+1)-ZSTAR(IET)) / DELTAZ
     &                * pas
!       ELSE
!         DZ(IPLOT) = 0.D0
!       ENDIF
!
        xp = xplot(iplot) + dx(iplot)
        yp = yplot(iplot) + dy(iplot)
        zp = zplot(iplot) + dz(iplot)
!
        shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                 -(y(i3)-y(i2))*(xp-x(i2))) * surdet(iel)
        shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                 -(y(i1)-y(i3))*(xp-x(i3))) * surdet(iel)
        shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                 -(y(i2)-y(i1))*(xp-x(i1))) * surdet(iel)
        shz(iplot) = (zp-zstar(iet)) / (zstar(iet+1)-zstar(iet))
!
        xplot(iplot) = xp
        yplot(iplot) = yp
        zplot(iplot) = zp
!
        IF(add) THEN
!         CONTINUOUS SETTING OF THE REACHED POSITION FOR IPLOT
!         AND THE NUMBER OF STEPS DONE ALREADY
          recvchar(iplot)%XP=xplot(iplot)
          recvchar(iplot)%YP=yplot(iplot)
          recvchar(iplot)%ZP=zplot(iplot)
          recvchar(iplot)%DX=dx(iplot)
          recvchar(iplot)%DY=dy(iplot)
          recvchar(iplot)%DZ=dz(iplot)
          recvchar(iplot)%INE=elt(iplot)
          recvchar(iplot)%ISP=isp
        ENDIF
!
!-----------------------------------------------------------------------
!       TEST: IS THE PATHLINE WENT OUT THE ORIGINAL ELEMENT
!-----------------------------------------------------------------------
!
50      CONTINUE
!
        iso = 0
        IF(shp(1,iplot).LT.epsilo)      iso=ibset(iso,2)
        IF(shp(2,iplot).LT.epsilo)      iso=ibset(iso,3)
        IF(shp(3,iplot).LT.epsilo)      iso=ibset(iso,4)
        IF(shz(iplot)  .LT.epsilo)      iso=ibset(iso,0)
        IF(shz(iplot)  .GT.1.d0-epsilo) iso=ibset(iso,1)
!
        IF(iso.NE.0) THEN
!
!-----------------------------------------------------------------------
!         HERE WE ARE OUT OF THE ELEMENT
!-----------------------------------------------------------------------
!
          isoh = iand(iso,28)
          isov = iand(iso, 3)
          iel = elt(iplot)
          iet = eta(iplot)
          xp = xplot(iplot)
          yp = yplot(iplot)
          zp = zplot(iplot)
!
          IF(isoh.NE.0) THEN
!
            IF(isoh.EQ.4) THEN
              ifa = 2
            ELSEIF(isoh.EQ.8) THEN
              ifa = 3
            ELSEIF(isoh.EQ.16) THEN
              ifa = 1
            ELSEIF(isoh.EQ.12) THEN
              ifa = 2
              IF(dx(iplot)*(y(ikle2(iel,3))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,3))-xp)) ifa = 3
            ELSEIF(isoh.EQ.24) THEN
              ifa = 3
              IF(dx(iplot)*(y(ikle2(iel,1))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,1))-xp)) ifa = 1
            ELSE
              ifa = 1
              IF(dx(iplot)*(y(ikle2(iel,2))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,2))-xp)) ifa = 2
            ENDIF
!
            IF(isov.GT.0) THEN
              IF(abs(dz(iplot)).GT.epsdz) THEN
                a1 = (zp-zstar(iet+isov-1)) / dz(iplot)
              ELSE
                a1 = 0.d0
              ENDIF
              i1 = ikle2(iel,ifa)
              i2 = ikle2(iel,isui(ifa))
              IF ((x(i2)-x(i1))*(yp-a1*dy(iplot)-y(i1)).GT.
     &            (y(i2)-y(i1))*(xp-a1*dx(iplot)-x(i1))) ifa=isov+3
            ENDIF
!
          ELSE
!
            ifa = isov + 3
!
          ENDIF
!
!         PERIODICITY (ALL VALUES ARE THE SAME SO FAR ON THE VERTICAL
!                      BUT WITH PERIODICITY WE MAY HAVE IET=NPLAN AND
!                      IT IS NOT POSSIBLE IN IFABOR).
!         IEL = IBOR(IEL,IFA,IET)
          iel = ibor(iel,ifa,1)
!         END PERIODICITY
!
          IF(ifa.LE.3) THEN
!
!-----------------------------------------------------------------------
!           HERE WE ARRIVE IN ANOTHER ELEMENT THROUGH A QUADRANGULAR FACE
!-----------------------------------------------------------------------
!
            IF(iel.GT.0) THEN
!
!-----------------------------------------------------------------------
!             RELOCALISING IN ADJACENT ELEMENT
!-----------------------------------------------------------------------
!
              i1 = ikle2(iel,1)
              i2 = ikle2(iel,2)
              i3 = ikle2(iel,3)
              elt(iplot) = iel
              shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
              shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
              shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
!
              GOTO 50
!
            ENDIF
!
!-----------------------------------------------------------------------
!           HERE WE PASS TO A NEIGHBOUR SUBDOMAIN AND COLLECT DATA
!-----------------------------------------------------------------------
!
            IF(iel.EQ.-2) THEN
              IF(.NOT.add) THEN
!               INTERFACE CROSSING
                CALL collect_char
     &            (ipid,iplot,elt(iplot),ifa,eta(iplot),0,isp,
     &             nsp,xplot(iplot),yplot(iplot),
     &             zplot(iplot),0.d0,
     &             dx(iplot),dy(iplot),dz(iplot),0.d0,
     &             ifapar,nchdim,nchara)
              ELSE
!               A LOST-AGAIN TRACEBACK DETECTED
!               PROCESSOR NUMBER
                recvchar(iplot)%NEPID=ifapar(ifa,elt(iplot))
                recvchar(iplot)%INE=ifapar(ifa+3,elt(iplot))
                recvchar(iplot)%KNE=eta(iplot)
              ENDIF
!             EXITING LOOP ON ISP
              EXIT
            ENDIF
!
!-----------------------------------------------------------------------
!           SPECIAL TREATMENT FOR SOLID OR LIQUID BOUNDARIES
!-----------------------------------------------------------------------
!
            dxp = dx(iplot)
            dyp = dy(iplot)
            i1  = ikle2(elt(iplot),ifa)
            i2  = ikle2(elt(iplot),isui(ifa))
            dx1 = x(i2) - x(i1)
            dy1 = y(i2) - y(i1)
!
            IF(iel.EQ.-1) THEN
!
!-----------------------------------------------------------------------
!             HERE SOLID BOUNDARY, VELOCITY IS PROJECTED ON THE BOUNDARY
!             AND WE GO ON
!-----------------------------------------------------------------------
!
!  TOMAWAC DIFFERENCES
!
!  STANDARD IMPLEMENTATION
!
!!            A1 = (DXP*DX1+DYP*DY1) / (DX1**2+DY1**2)
!!            DX(IPLOT) = A1 * DX1
!!            DY(IPLOT) = A1 * DY1
!
!!            A1=((XP-X(I1))*DX1+(YP-Y(I1))*DY1)/(DX1**2+DY1**2)
!!            SHP(          IFA  ,IPLOT) = 1.D0 - A1
!!            SHP(     ISUI(IFA) ,IPLOT) = A1
!!            SHP(ISUI(ISUI(IFA)),IPLOT) = 0.D0
!!            XPLOT(IPLOT) = X(I1) + A1 * DX1
!!            YPLOT(IPLOT) = Y(I1) + A1 * DY1
!!            IF(ADD) THEN
!!              RECVCHAR(IPLOT)%XP=XPLOT(IPLOT)
!!              RECVCHAR(IPLOT)%YP=YPLOT(IPLOT)
!!              RECVCHAR(IPLOT)%ZP=ZPLOT(IPLOT)
!!              RECVCHAR(IPLOT)%DX=DX(IPLOT)
!!              RECVCHAR(IPLOT)%DY=DY(IPLOT)
!!              RECVCHAR(IPLOT)%DZ=DZ(IPLOT)
!!            ENDIF
!
!!            GOTO 50
!
!  TOMAWAC IMPLEMENTATION
!
              shp(1,iplot) = 0.d0
              shp(2,iplot) = 0.d0
              shp(3,iplot) = 0.d0
              elt(iplot) = - sens * elt(iplot)
              EXIT  ! LOOP ON ISP
!
            ELSEIF(iel.EQ.0) THEN
!
!-----------------------------------------------------------------------
!           HERE WE HAVE A LIQUID BOUNDARY, THE CHARACTERISTIC IS STOPPED
!-----------------------------------------------------------------------
!
            denom = dxp*dy1-dyp*dx1
            IF(abs(denom).GT.1.d-12) THEN
              a1 = (dxp*(yp-y(i1))-dyp*(xp-x(i1))) / denom
            ELSE
              a1 = 0.d0
            ENDIF
            IF(a1.GT.1.d0) a1 = 1.d0
            IF(a1.LT.0.d0) a1 = 0.d0
            shp(          ifa  ,iplot) = 1.d0 - a1
            shp(     isui(ifa) ,iplot) = a1
            shp(isui(isui(ifa)),iplot) = 0.d0
            xplot(iplot) = x(i1) + a1 * dx1
            yplot(iplot) = y(i1) + a1 * dy1
            IF(abs(dxp).GT.abs(dyp)) THEN
              a1 = (xp-xplot(iplot))/dxp
            ELSE
              a1 = (yp-yplot(iplot))/dyp
            ENDIF
            zplot(iplot) = zp - a1*dz(iplot)
            shz(iplot) = (zplot(iplot)-zstar(iet))
     &                 / (zstar(iet+1)-zstar(iet))
!           THIS IS A MARKER FOR PARTICLES EXITING A DOMAIN
!           SENS=-1 FOR BACKWARD CHARACTERISTICS
            elt(iplot) = - sens * elt(iplot)
!           EXITING LOOP ON ISP
            EXIT
!
            ELSE
!
              WRITE(lu,*) 'UNEXPECTED CASE IN SCHAR41'
              WRITE(lu,*) 'IEL=',iel
              CALL plante(1)
              stop
!
            ENDIF
!
          ELSE
!
!-----------------------------------------------------------------------
!  CASE IFA = 4 OR 5
!  HERE WE EXIT THROUGH TOP OR BOTTOM OF THE PRISM
!-----------------------------------------------------------------------
!
            ifa = ifa - 4
!           HENCE IFA NOW EQUALS 0 OR 1
!
            IF(iel.EQ.1) THEN
!
!-----------------------------------------------------------------------
!           NO NEED TO RECOMPUTE THE VELOCITIES,
!           RELOCALISING IN NEW ELEMENT
!-----------------------------------------------------------------------
!
              eta(iplot) = iet + ifa + ifa - 1
!! PERIODICITY (THIS CAN NEVER HAPPEN WITHOUT PERIODICITY)
              IF(eta(iplot).EQ.nplan+1) THEN
                eta(iplot)=1
                zp=zp-zstar(nplan+1)
                zplot(iplot)=zp
              ENDIF
              IF(eta(iplot).EQ.0) THEN
                eta(iplot) = nplan
                zp=zp+zstar(nplan+1)
                zplot(iplot)=zp
              ENDIF
!! END OF PERIODICITY
              shz(iplot) = (zp-zstar(eta(iplot)))
     &                   / (zstar(eta(iplot)+1)-zstar(eta(iplot)))
!
              IF(add) THEN
                recvchar(iplot)%KNE=eta(iplot)
                recvchar(iplot)%ZP=zp
              ENDIF
!
              iso = isoh
!
              IF(shz(iplot).LT.     epsilo) iso=ibset(iso,0)
              IF(shz(iplot).GT.1.d0-epsilo) iso=ibset(iso,1)
!
              GOTO 50
!
            ENDIF
!
            IF(iel.EQ.-1) THEN
!
!-----------------------------------------------------------------------
!  LA, ON SAIT QUE LA FACE DE SORTIE EST UNE FRONTIERE SOLIDE
!  ON PROJETTE LE RELIQUAT SUR LA FRONTIERE PUIS ON SE RELOCALISE
!-----------------------------------------------------------------------
!
              zplot(iplot) = zstar(iet+ifa)
              dz(iplot) = 0.d0
              shz(iplot) = ifa
!
              iso = isoh
              IF(isoh.NE.0) GOTO 50
!
            ELSE
!
!-----------------------------------------------------------------------
!  LA, ON SAIT QUE LA FACE DE SORTIE EST UNE FRONTIERE LIQUIDE (CAS 0)
!      ON ARRETE ALORS LA REMONTEE DES CARACTERISTIQUES (SIGNE DE ELT)
!  OU, QUE L'ON VIENT DE TRAVERSER UN PLAN AVEC RECALCUL DES VITESSES
!  DEMANDE (CAS 2)
!-----------------------------------------------------------------------
!
              IF(abs(dz(iplot)).GT.epsdz) THEN
                a1 = (zp-zstar(iet+ifa)) / dz(iplot)
              ELSE
                a1 = 0.d0
              ENDIF
              xp = xp - a1*dx(iplot)
              yp = yp - a1*dy(iplot)
              zp = zstar(iet+ifa)
              iele = elt(iplot)
              i1 = ikle2(iele,1)
              i2 = ikle2(iele,2)
              i3 = ikle2(iele,3)
!
              shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iele)
              shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iele)
              shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iele)
!
              IF(iel.EQ.2) THEN
!
!-----------------------------------------------------------------------
!  LA, ON SAIT QUE LA FACE DE SORTIE SE SITUE SUR LE PLAN OU ON DEMANDE
!  UN RECALCUL DES VITESSES
!-----------------------------------------------------------------------
!
!               IF IFA = 1 EXIT THROUGH THE TOP
!               IF IFA = 0 EXIT THROUGH THE BOTTOM
!               THEN NEW IET IS  IET+1 IF IFA = 1
!                            AND IET-1 IF IFA = 0
!               THIS IS SUMMARISED BY IET=IET+2*IFA-1
!
!               RECOMPUTED VELOCITIES MUST BE TAKEN AT IET2=IET+IFA
!               I.E. BOTTOM IF EXIT THROUGH THE BOTTOM
!               AND TOP IF EXIT THROUGH THE TOP
!
                iet2 = iet + ifa
                iet  = iet + ifa + ifa - 1
!
                pas2 = pas2 * a1
!
                dx(iplot) = ( u(i1,iet2)*shp(1,iplot)
     &                      + u(i2,iet2)*shp(2,iplot)
     &                      + u(i3,iet2)*shp(3,iplot) ) * pas2
!
                dy(iplot) = ( v(i1,iet2)*shp(1,iplot)
     &                      + v(i2,iet2)*shp(2,iplot)
     &                      + v(i3,iet2)*shp(3,iplot) ) * pas2
!
!               DELTAZ =  (Z(I1,IET+1)-Z(I1,IET))*SHP(1,IPLOT)
!    &                  + (Z(I2,IET+1)-Z(I2,IET))*SHP(2,IPLOT)
!    &                  + (Z(I3,IET+1)-Z(I3,IET))*SHP(3,IPLOT)
!
!               IF(DELTAZ.GT.EPSDZ) THEN
                  dz(iplot) = ( w(i1,iet2)*shp(1,iplot)
     &                        + w(i2,iet2)*shp(2,iplot)
     &                        + w(i3,iet2)*shp(3,iplot) ) * pas2
!    &                        * (ZSTAR(IET+1)-ZSTAR(IET)) / DELTAZ
!               ELSE
!                 DZ(IPLOT) = 0.D0
!               ENDIF
!
                xp = xp + dx(iplot)
                yp = yp + dy(iplot)
                zp = zp + dz(iplot)
!
                shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2))) * surdet(iele)
                shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3))) * surdet(iele)
                shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1))) * surdet(iele)
                shz(iplot)=(zp-zstar(iet))/(zstar(iet+1)-zstar(iet))
!
                xplot(iplot) = xp
                yplot(iplot) = yp
                zplot(iplot) = zp
                eta(iplot) = iet
!
                iso = 0
!
                IF(shp(1,iplot).LT.epsilo) iso=ibset(iso,2)
                IF(shp(2,iplot).LT.epsilo) iso=ibset(iso,3)
                IF(shp(3,iplot).LT.epsilo) iso=ibset(iso,4)
!
                IF(shz(iplot).LT.     epsilo) iso=ibset(iso,0)
                IF(shz(iplot).GT.1.d0-epsilo) iso=ibset(iso,1)
!
                IF(add) THEN
                  recvchar(iplot)%XP=xplot(iplot)
                  recvchar(iplot)%YP=yplot(iplot)
                  recvchar(iplot)%ZP=zplot(iplot)
                  recvchar(iplot)%DX=dx(iplot)
                  recvchar(iplot)%DY=dy(iplot)
                  recvchar(iplot)%DZ=dz(iplot)
                ENDIF
!
                GOTO 50
!
              ENDIF
!
              xplot(iplot) = xp
              yplot(iplot) = yp
              zplot(iplot) = zp
              shz(iplot) = ifa
!             THIS IS A MARKER FOR PARTICLES EXITING A DOMAIN
!             SENS=-1 FOR BACKWARD CHARACTERISTICS
              elt(iplot) = - sens * elt(iplot)
!             EXITING LOOP ON ISP
              EXIT
!
            ENDIF
!
          ENDIF
!
!       IF(ISO.NE.0)
        ENDIF
!
        ENDDO
      ENDDO
!
      IF(.NOT.add) THEN
        IF(ncsize.GT.1) nspmax=p_max(nspmax)
        WRITE(lu,*) 'NUMBER OF SUB-ITERATIONS :',nspmax
      ENDIF
!
!     RESTORING ORIGINAL ETA1
!
      eta1(nplan)=nplan+1
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE schar41_per
!                       *************************
                        SUBROUTINE schar41_per_4d
!                       *************************
!
     &( u , v , w , f , dt , nrk , x , y , zstar , freq ,ikle2 ,ibor ,
     &  xplot , yplot , zplot  , fplot,dx , dy , dz , df,
     &  shp   , shz   , shf    , elt    , eta ,
     &  fre   , nplot , npoin2 , nelem2 , nelmax2,nplan , nf,surdet ,
     &  sens  , ifapar, nchdim ,nchara,add)
!
!***********************************************************************
! BIEF VERSION 7.3
!
!brief    LIKE SCHAR41 BUT WITH PERIODICITY ON THE VERTICAL
!         DIFFERENCES MAEKED WITH "PERIODICITY"
!
!
!history  J-M HERVOUET (LNHE)
!+        12/07/2012
!+        V6P3
!+      First version (differences with schar41 taken in Tomawac)
!
!***********************************************************************
!
!  FONCTION :
!
!     REMONTEE OU DESCENTE
!     DES COURBES CARACTERISTIQUES
!     SUR DES PRISMES DE TELEMAC-3D
!     DANS L'INTERVALLE DE TEMPS DT
!     AVEC UNE DISCRETISATION ELEMENTS FINIS
!
!
!  DISCRETISATION :
!
!     LE DOMAINE EST APPROCHE PAR UNE DISCRETISATION ELEMENTS FINIS
!     UNE APPROXIMATION LOCALE EST DEFINIE POUR LE VECTEUR VITESSE :
!     LA VALEUR EN UN POINT D'UN ELEMENT NE DEPEND QUE DES VALEURS
!     AUX NOEUDS DE CET ELEMENT
!
!
!  RESTRICTIONS ET HYPOTHESES :
!
!     LE CHAMP CONVECTEUR U EST SUPPOSE INDEPENDANT DU TEMPS
!
!-----------------------------------------------------------------------
!                             ARGUMENTS
! .________________.____.______________________________________________.
! |      NOM       |MODE|                   ROLE                       |
! |________________|____|______________________________________________|
! |    U,V,W       | -->| COMPONENTS OF ADVECTION VELOCITY
! |                |    | BUT W IS W* x DELTAZ (STEMS FROM TRIDW2)
! |    DT          | -->| PAS DE TEMPS.                                |
! |    NRK         | -->| NOMBRE DE SOUS-PAS DE RUNGE-KUTTA.           |
! |    X,Y,ZSTAR   | -->| COORDONNEES DES POINTS DU MAILLAGE.          | |
! |    IKLE2       | -->| TRANSITION ENTRE LES NUMEROTATIONS LOCALE    |
! |                |    | ET GLOBALE DU MAILLAGE 2D.                   |
! |    IBOR        | -->| NUMEROS 2D DES ELEMENTS AYANT UNE FACE COMMUNE
! |                |    | AVEC L'ELEMENT .  SI IBOR<0 OU NUL           |
! |                |    | ON A UNE FACE LIQUIDE,SOLIDE,OU PERIODIQUE   |
! |  X..,Y..,ZPLOT |<-->| POSITIONS SUCCESSIVES DES DERIVANTS.         |
! |    DX,DY,DZ    | -- | STOCKAGE DES SOUS-PAS .                      |
! |    SHP         |<-->| COORDONNEES BARYCENTRIQUES 2D AU PIED DES    |
! |                |    | COURBES CARACTERISTIQUES.                    |
! |    SHZ         |<-->| COORDONNEES BARYCENTRIQUES SUIVANT Z DES     |
! |                |    | NOEUDS DANS LEURS ETAGES "ETA" ASSOCIES.     |
! |    ELT         |<-->| NUMEROS DES ELEMENTS 2D CHOISIS POUR CHAQUE  |
! |                |    | NOEUD.                                       |
! |    ETA         |<-->| NUMEROS DES ETAGES CHOISIS POUR CHAQUE NOEUD.|
! |    NPLOT       | -->| NOMBRE DE DERIVANTS.                         |
! |    NPOIN2      | -->| NOMBRE DE POINTS DU MAILLAGE 2D.             |
! |    NELEM2      | -->| NOMBRE D'ELEMENTS DU MAILLAGE 2D.            |
! |    NPLAN       | -->| NOMBRE DE PLANS.                             |
! |    SURDET      | -->| VARIABLE UTILISEE PAR LA TRANSFORMEE ISOPARAM.
! |    SENS        | -->| DESCENTE OU REMONTEE DES CARACTERISTIQUES.   |
! |    ISO         | -- | STOCKAGE BINAIRE DE LA FACE DE SORTIE.       |
! |________________|____|______________________________________________|
!  MODE: -->(DONNEE NON MODIFIEE),<--(RESULTAT),<-->(DONNEE MODIFIEE)
!-----------------------------------------------------------------------
!     - APPELE PAR : CARACT
!     - PROGRAMMES APPELES : NEANT
!
!***********************************************************************
!
      USE bief
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)    :: SENS,NPLAN,NCHDIM,NF,NELMAX2
      INTEGER         , INTENT(IN)    :: NPOIN2,NELEM2,NPLOT,NRK
      INTEGER         , INTENT(IN)    :: IKLE2(nelmax2,3)
      INTEGER         , INTENT(INOUT) :: ELT(nplot),NCHARA
      DOUBLE PRECISION, INTENT(IN)    :: U(npoin2,nplan,nf)
      DOUBLE PRECISION, INTENT(IN)    :: V(npoin2,nplan,nf)
      DOUBLE PRECISION, INTENT(IN)    :: W(npoin2,nplan,nf)
      DOUBLE PRECISION, INTENT(IN)    :: F(npoin2,nplan,nf)
      DOUBLE PRECISION, INTENT(IN)    :: SURDET(nelem2)
      DOUBLE PRECISION, INTENT(INOUT) :: XPLOT(nplot),YPLOT(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: ZPLOT(nplot),FPLOT(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: SHP(3,nplot),SHZ(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: SHF(nplot)
      DOUBLE PRECISION, INTENT(IN)    :: X(npoin2),Y(npoin2),DT
!                                              PERIODICITY
      DOUBLE PRECISION, INTENT(IN)    :: ZSTAR(nplan+1)
      DOUBLE PRECISION, INTENT(IN)    :: FREQ(nf)
      DOUBLE PRECISION, INTENT(INOUT) :: DX(nplot),DY(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: DZ(nplot),DF(nplot)
      INTEGER         , INTENT(IN)    :: IBOR(nelmax2,5,nplan-1)
      INTEGER         , INTENT(INOUT) :: ETA(nplot),FRE(nplot)
      INTEGER         , INTENT(IN)    :: IFAPAR(6,*)
      LOGICAL, INTENT(IN)             :: ADD
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER ISO,ISPDONE,NSP,NSPMAX
      INTEGER IPLOT,ISP,I1,I2,I3,IEL,IET,ISOH,ISOV,ISOF,ISOT
      INTEGER :: IETP1,IFA,ISUI(3),IFR
!
      DOUBLE PRECISION PAS,A1,A2,DX1,DY1,DXP,DYP,DZP,XP,YP,ZP,FP
      DOUBLE PRECISION PAS2,DFP,DENOM
!
      INTRINSIC abs , int , max , sqrt
!
      parameter( isui = (/ 2 , 3 , 1 /) )
      DOUBLE PRECISION, PARAMETER :: EPSILO = -1.d-6
!
      nspmax=1
      eta1(nplan)=1
!
!-----------------------------------------------------------------------
!     FOR EVERY POINT
!-----------------------------------------------------------------------
!
      DO iplot = 1 , nplot
!
        IF(add) THEN
!
        xplot(iplot)   = recvchar(iplot)%XP
        yplot(iplot)   = recvchar(iplot)%YP
        zplot(iplot)   = recvchar(iplot)%ZP
        fplot(iplot)   = recvchar(iplot)%FP
        dx(iplot)      = recvchar(iplot)%DX
        dy(iplot)      = recvchar(iplot)%DY
        dz(iplot)      = recvchar(iplot)%DZ
        df(iplot)      = recvchar(iplot)%DF
        elt(iplot)     = recvchar(iplot)%INE
        eta(iplot)     = recvchar(iplot)%KNE
        fre(iplot)     = recvchar(iplot)%IFR
        nsp            = recvchar(iplot)%NSP ! R-K STEPS TO BE FULLFILLED
        ispdone        = recvchar(iplot)%ISP ! R-K STEPS ALREADY DONE
        iel = elt(iplot)
        iet = eta(iplot)
        ifr = fre(iplot)
        xp  = xplot(iplot)
        yp  = yplot(iplot)
        zp  = zplot(iplot)
        fp  = fplot(iplot)
        i1 = ikle2(iel,1)
        i2 = ikle2(iel,2)
        i3 = ikle2(iel,3)
        shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                 -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
        shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                 -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
        shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                 -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
        shz(iplot) = (zp-zstar(iet)) / (zstar(iet+1)-zstar(iet))
        shf(iplot) = (fp-freq(ifr)) / (freq(ifr+1)-freq(ifr))
!       ASSUME ALL ARE LOCALISED, IT WILL BE SET OTHERWISE IF LOST-AGAIN
        recvchar(iplot)%NEPID=-1
!
        ELSE
!
          iel = elt(iplot)
!         POINTS WITH IEL=0 ARE TREATED SO THAT THE FINAL
!         INTERPOLATION GIVES 0.,
!         AND WE SKIP TO NEXT POINT IPLOT (CYCLE)
          IF(iel.EQ.0) THEN
            elt(iplot)=1
            eta(iplot)=1
            fre(iplot)=1
            shp(1,iplot)=0.d0
            shp(2,iplot)=0.d0
            shp(3,iplot)=0.d0
            shz(iplot)=0.d0
            shf(iplot)=0.d0
            cycle
          ENDIF
          iet = eta(iplot)
          ifr = fre(iplot)
          i1 = ikle2(iel,1)
          i2 = ikle2(iel,2)
          i3 = ikle2(iel,3)
!         HERE IET+1 IS ALWAYS < NPLAN+1 (SEE GTSH41)
          dxp =(1.d0-shf(iplot))*
     &           ((u(i1,iet  ,ifr)*shp(1,iplot)
     &           + u(i2,iet  ,ifr)*shp(2,iplot)
     &           + u(i3,iet  ,ifr)*shp(3,iplot))*(1.d0-shz(iplot))
     &           +(u(i1,iet+1,ifr)*shp(1,iplot)
     &           + u(i2,iet+1,ifr)*shp(2,iplot)
     &           + u(i3,iet+1,ifr)*shp(3,iplot))*shz(iplot))
     &         + shf(iplot)*
     &           ((u(i1,iet  ,ifr+1)*shp(1,iplot)
     &           + u(i2,iet  ,ifr+1)*shp(2,iplot)
     &           + u(i3,iet  ,ifr+1)*shp(3,iplot))*(1.d0-shz(iplot))
     &           +(u(i1,iet+1,ifr+1)*shp(1,iplot)
     &           + u(i2,iet+1,ifr+1)*shp(2,iplot)
     &           + u(i3,iet+1,ifr+1)*shp(3,iplot))*shz(iplot))
!
          dyp =(1.d0-shf(iplot))*
     &           ((v(i1,iet  ,ifr)*shp(1,iplot)
     &           + v(i2,iet  ,ifr)*shp(2,iplot)
     &           + v(i3,iet  ,ifr)*shp(3,iplot))*(1.d0-shz(iplot))
     &           +(v(i1,iet+1,ifr)*shp(1,iplot)
     &           + v(i2,iet+1,ifr)*shp(2,iplot)
     &           + v(i3,iet+1,ifr)*shp(3,iplot))*shz(iplot))
     &         + shf(iplot)*
     &           ((v(i1,iet  ,ifr+1)*shp(1,iplot)
     &           + v(i2,iet  ,ifr+1)*shp(2,iplot)
     &           + v(i3,iet  ,ifr+1)*shp(3,iplot))*(1.d0-shz(iplot))
     &           +(v(i1,iet+1,ifr+1)*shp(1,iplot)
     &           + v(i2,iet+1,ifr+1)*shp(2,iplot)
     &           + v(i3,iet+1,ifr+1)*shp(3,iplot))*shz(iplot))
!
          dzp =(1.d0-shf(iplot))*
     &           ((w(i1,iet  ,ifr)*shp(1,iplot)
     &           + w(i2,iet  ,ifr)*shp(2,iplot)
     &           + w(i3,iet  ,ifr)*shp(3,iplot))*(1.d0-shz(iplot))
     &           +(w(i1,iet+1,ifr)*shp(1,iplot)
     &           + w(i2,iet+1,ifr)*shp(2,iplot)
     &           + w(i3,iet+1,ifr)*shp(3,iplot))*shz(iplot))
     &         + shf(iplot)*
     &           ((w(i1,iet  ,ifr+1)*shp(1,iplot)
     &           + w(i2,iet  ,ifr+1)*shp(2,iplot)
     &           + w(i3,iet  ,ifr+1)*shp(3,iplot))*(1.d0-shz(iplot))
     &           +(w(i1,iet+1,ifr+1)*shp(1,iplot)
     &           + w(i2,iet+1,ifr+1)*shp(2,iplot)
     &           + w(i3,iet+1,ifr+1)*shp(3,iplot))*shz(iplot))
!
          dfp =(1.d0-shf(iplot))*
     &           ((f(i1,iet  ,ifr)*shp(1,iplot)
     &           + f(i2,iet  ,ifr)*shp(2,iplot)
     &           + f(i3,iet  ,ifr)*shp(3,iplot))*(1.d0-shz(iplot))
     &           +(f(i1,iet+1,ifr)*shp(1,iplot)
     &           + f(i2,iet+1,ifr)*shp(2,iplot)
     &           + f(i3,iet+1,ifr)*shp(3,iplot))*shz(iplot))
     &         + shf(iplot)*
     &           ((f(i1,iet  ,ifr+1)*shp(1,iplot)
     &           + f(i2,iet  ,ifr+1)*shp(2,iplot)
     &           + f(i3,iet  ,ifr+1)*shp(3,iplot))*(1.d0-shz(iplot))
     &           +(f(i1,iet+1,ifr+1)*shp(1,iplot)
     &           + f(i2,iet+1,ifr+1)*shp(2,iplot)
     &           + f(i3,iet+1,ifr+1)*shp(3,iplot))*shz(iplot))
!
          nsp=max(1,int(nrk*dt*sqrt((dxp**2+dyp**2)*surdet(iel))))
          nsp=max(nsp,int(nrk*dt*abs(dzp/(zstar(iet+1)-zstar(iet)))))
          nsp=max(nsp,int(nrk*dt*abs(dfp/(freq(ifr)-freq(ifr+1)))))
          nspmax=max(nsp,nspmax)
!!        END PERIODICITY
!
          ispdone=1
!
        ENDIF
!
        pas = sens * dt / nsp
!
!       LOOP ON RUNGE-KUTTA SUB-STEPS
!
!       COMPILER MUST DO NOTHING IF ISPDONE>NSP
!       IN MODE "ADD", ISP = ISPDONE HAS NOT BEEN FULLY DONE
!       IT IS RESTARTED HERE
!
      DO isp = ispdone,nsp
!
!-----------------------------------------------------------------------
!       LOCALISING THE ARRIVAL POINT
!-----------------------------------------------------------------------
!
        pas2 = pas
!
!       IN MODE "ADD" ITERATIONS ALREADY DONE ARE SKIPPED AND
!                     CHARACTERISTICS GONE IN ANOTHER SUB-DOMAIN SKIPPED
!
        IF(add) THEN
          IF(isp.EQ.ispdone) GO TO 50
          IF(recvchar(iplot)%NEPID.NE.-1) cycle
        ENDIF
!
        iel = elt(iplot)
        iet = eta(iplot)
        ifr = fre(iplot)
!!      PERIODICITY IETP1 REPLACES IET+1 (BUT NOT ALWAYS)
        ietp1=eta1(iet)
        i1 = ikle2(iel,1)
        i2 = ikle2(iel,2)
        i3 = ikle2(iel,3)
!
        dx(iplot) = ( (1.d0-shf(iplot))*
     &      ((u(i1,iet  ,ifr)*shp(1,iplot)
     &      + u(i2,iet  ,ifr)*shp(2,iplot)
     &      + u(i3,iet  ,ifr)*shp(3,iplot))*(1.d0-shz(iplot))
     &      +(u(i1,ietp1,ifr)*shp(1,iplot)
     &      + u(i2,ietp1,ifr)*shp(2,iplot)
     &      + u(i3,ietp1,ifr)*shp(3,iplot))*shz(iplot))
     &        + shf(iplot)*
     &      ((u(i1,iet  ,ifr+1)*shp(1,iplot)
     &      + u(i2,iet  ,ifr+1)*shp(2,iplot)
     &      + u(i3,iet  ,ifr+1)*shp(3,iplot))*(1.d0-shz(iplot))
     &      +(u(i1,ietp1,ifr+1)*shp(1,iplot)
     &      + u(i2,ietp1,ifr+1)*shp(2,iplot)
     &      + u(i3,ietp1,ifr+1)*shp(3,iplot))*shz(iplot)))*pas
!
        dy(iplot) = ( (1.d0-shf(iplot))*
     &      ((v(i1,iet  ,ifr)*shp(1,iplot)
     &      + v(i2,iet  ,ifr)*shp(2,iplot)
     &      + v(i3,iet  ,ifr)*shp(3,iplot))*(1.d0-shz(iplot))
     &      +(v(i1,ietp1,ifr)*shp(1,iplot)
     &      + v(i2,ietp1,ifr)*shp(2,iplot)
     &      + v(i3,ietp1,ifr)*shp(3,iplot))*shz(iplot))
     &        + shf(iplot)*
     &      ((v(i1,iet  ,ifr+1)*shp(1,iplot)
     &      + v(i2,iet  ,ifr+1)*shp(2,iplot)
     &      + v(i3,iet  ,ifr+1)*shp(3,iplot))*(1.d0-shz(iplot))
     &      +(v(i1,ietp1,ifr+1)*shp(1,iplot)
     &      + v(i2,ietp1,ifr+1)*shp(2,iplot)
     &      + v(i3,ietp1,ifr+1)*shp(3,iplot))*shz(iplot)))*pas
!
        dz(iplot) = ( (1.d0-shf(iplot))*
     &      ((w(i1,iet  ,ifr)*shp(1,iplot)
     &      + w(i2,iet  ,ifr)*shp(2,iplot)
     &      + w(i3,iet  ,ifr)*shp(3,iplot))*(1.d0-shz(iplot))
     &      +(w(i1,ietp1,ifr)*shp(1,iplot)
     &      + w(i2,ietp1,ifr)*shp(2,iplot)
     &      + w(i3,ietp1,ifr)*shp(3,iplot))*shz(iplot))
     &        + shf(iplot)*
     &      ((w(i1,iet  ,ifr+1)*shp(1,iplot)
     &      + w(i2,iet  ,ifr+1)*shp(2,iplot)
     &      + w(i3,iet  ,ifr+1)*shp(3,iplot))*(1.d0-shz(iplot))
     &      +(w(i1,ietp1,ifr+1)*shp(1,iplot)
     &      + w(i2,ietp1,ifr+1)*shp(2,iplot)
     &      + w(i3,ietp1,ifr+1)*shp(3,iplot))*shz(iplot)))*pas
!
        df(iplot) = ( (1.d0-shf(iplot))*
     &      ((f(i1,iet  ,ifr)*shp(1,iplot)
     &      + f(i2,iet  ,ifr)*shp(2,iplot)
     &      + f(i3,iet  ,ifr)*shp(3,iplot))*(1.d0-shz(iplot))
     &      +(f(i1,ietp1,ifr)*shp(1,iplot)
     &      + f(i2,ietp1,ifr)*shp(2,iplot)
     &      + f(i3,ietp1,ifr)*shp(3,iplot))*shz(iplot))
     &        + shf(iplot)*
     &      ((f(i1,iet  ,ifr+1)*shp(1,iplot)
     &      + f(i2,iet  ,ifr+1)*shp(2,iplot)
     &      + f(i3,iet  ,ifr+1)*shp(3,iplot))*(1.d0-shz(iplot))
     &      +(f(i1,ietp1,ifr+1)*shp(1,iplot)
     &      + f(i2,ietp1,ifr+1)*shp(2,iplot)
     &      + f(i3,ietp1,ifr+1)*shp(3,iplot))*shz(iplot)) )*pas
!
        xp = xplot(iplot) + dx(iplot)
        yp = yplot(iplot) + dy(iplot)
        zp = zplot(iplot) + dz(iplot)
        fp = fplot(iplot) + df(iplot)
!
        shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                 -(y(i3)-y(i2))*(xp-x(i2))) * surdet(iel)
        shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                 -(y(i1)-y(i3))*(xp-x(i3))) * surdet(iel)
        shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                 -(y(i2)-y(i1))*(xp-x(i1))) * surdet(iel)
        shz(iplot) = (zp-zstar(iet)) / (zstar(iet+1)-zstar(iet))
        shf(iplot) = (fp-freq(ifr)) / (freq(ifr+1)-freq(ifr))
!
        xplot(iplot) = xp
        yplot(iplot) = yp
        zplot(iplot) = zp
        fplot(iplot) = fp
!
        IF(add) THEN
!         CONTINUOUS SETTING OF THE REACHED POSITION FOR IPLOT
!         AND THE NUMBER OF STEPS DONE ALREADY
          recvchar(iplot)%XP=xplot(iplot)
          recvchar(iplot)%YP=yplot(iplot)
          recvchar(iplot)%ZP=zplot(iplot)
          recvchar(iplot)%FP=fplot(iplot)
          recvchar(iplot)%DX=dx(iplot)
          recvchar(iplot)%DY=dy(iplot)
          recvchar(iplot)%DZ=dz(iplot)
          recvchar(iplot)%DF=df(iplot)
          recvchar(iplot)%INE=elt(iplot)
          recvchar(iplot)%ISP=isp
        ENDIF
!
!-----------------------------------------------------------------------
!       TEST: IF THE PATHLINE WENT OUT THE ORIGINAL ELEMENT
!-----------------------------------------------------------------------
!
50      CONTINUE
!
        iso = 0
        IF(shp(1,iplot).LT.epsilo)    iso=ibset(iso,4)
        IF(shp(2,iplot).LT.epsilo)    iso=ibset(iso,5)
        IF(shp(3,iplot).LT.epsilo)    iso=ibset(iso,6)
        IF(shz(iplot).LT.epsilo)      iso=ibset(iso,0)
        IF(shz(iplot).GT.1.d0-epsilo) iso=ibset(iso,1)
        IF(shf(iplot).LT.epsilo)      iso=ibset(iso,2)
        IF(shf(iplot).GT.1.d0-epsilo) iso=ibset(iso,3)
!
        IF(iso.NE.0) THEN
!
!-----------------------------------------------------------------------
!         HERE: WE LEFT THE ELEMENT
!-----------------------------------------------------------------------
!
          isot = iand(iso, 3)
          isof = iand(iso,12)/4
          isov = iand(iso,15)
          isoh = iand(iso,112)
          iel  = elt(iplot)
          iet  = eta(iplot)
          ifr  = fre(iplot)
          xp   = xplot(iplot)
          yp   = yplot(iplot)
          zp   = zplot(iplot)
          fp   = fplot(iplot)
!
          IF(isoh.NE.0) THEN
!
            IF(isoh.EQ.16) THEN
              ifa = 2
            ELSEIF(isoh.EQ.32) THEN
              ifa = 3
            ELSEIF(isoh.EQ.64) THEN
              ifa = 1
            ELSEIF(isoh.EQ.48) THEN
              ifa = 2
              IF(dx(iplot)*(y(ikle2(iel,3))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,3))-xp)) ifa = 3
            ELSEIF(isoh.EQ.96) THEN
              ifa = 3
              IF(dx(iplot)*(y(ikle2(iel,1))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,1))-xp)) ifa = 1
            ELSE
!             CASE ISOH=80
              ifa = 1
              IF(dx(iplot)*(y(ikle2(iel,2))-yp).LT.
     &           dy(iplot)*(x(ikle2(iel,2))-xp)) ifa = 2
            ENDIF
!
            IF(isov.GT.0) THEN
              i1 = ikle2(iel,ifa)
              i2 = ikle2(iel,isui(ifa))
              IF(isof.GT.0) THEN
                IF(isot.GT.0) THEN
                  a1=(fp- freq(ifr+isof-1))/df(iplot)
                  a2=(zp-zstar(iet+isot-1))/dz(iplot)
                  IF(a1.LT.a2) THEN
                    IF((x(i2)-x(i1))*(yp-a1*dy(iplot)-y(i1)).GT.
     &                 (y(i2)-y(i1))*(xp-a1*dx(iplot)-x(i1))) THEN
                      ifa=isof+5
                    ENDIF
                  ELSE
                    IF((x(i2)-x(i1))*(yp-a2*dy(iplot)-y(i1)).GT.
     &             (y(i2)-y(i1))*(xp-a2*dx(iplot)-x(i1))) THEN
                      ifa=isot+3
                    ENDIF
                  ENDIF
                ELSE
                  a1 = (fp-freq(ifr+isof-1)) / df(iplot)
                  IF((x(i2)-x(i1))*(yp-a1*dy(iplot)-y(i1)).GT.
     &               (y(i2)-y(i1))*(xp-a1*dx(iplot)-x(i1))) THEN
                    ifa=isof+5
                  ENDIF
                ENDIF
              ELSE
                a1 = (zp-zstar(iet+isot-1)) / dz(iplot)
                IF((x(i2)-x(i1))*(yp-a1*dy(iplot)-y(i1)).GT.
     &             (y(i2)-y(i1))*(xp-a1*dx(iplot)-x(i1))) THEN
                  ifa=isot+3
                ENDIF
              ENDIF
            ENDIF
!
          ELSEIF(isot.GT.0) THEN
            ifa = isot + 3
            IF(isof.GT.0) THEN
              a1=(fp-freq(ifr+isof-1))/df(iplot)
              a2=(zp-zstar(iet+isot-1))/dz(iplot)
              IF(a1.LT.a2) ifa = isof + 5
            ENDIF
          ELSE
            ifa = isof + 5
          ENDIF
!
          IF(ifa.LE.3) THEN
!
!         IEL = IBOR(IEL,IFA,IET)
          iel = ibor(iel,ifa,1)
!
!-----------------------------------------------------------------------
!  HERE: THE EXIT FACE OF THE PRISM IS A RECTANGULAR FACE
!  =================================================================
!-----------------------------------------------------------------------
!
            IF(iel.GT.0) THEN
!
!-----------------------------------------------------------------------
!  HERE: THE EXIT FACE IS AN INTERIOR FACE
!  MOVES TO THE ADJACENT ELEMENT
!-----------------------------------------------------------------------
!
              i1 = ikle2(iel,1)
              i2 = ikle2(iel,2)
              i3 = ikle2(iel,3)
!
              elt(iplot) = iel
              shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
              shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
              shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
!
              GOTO 50
!
            ENDIF
!
!-----------------------------------------------------------------------
!           HERE WE PASS TO A NEIGHBOUR SUBDOMAIN AND COLLECT DATA
!-----------------------------------------------------------------------
!
            IF(iel.EQ.-2) THEN
              IF(.NOT.add) THEN
!               INTERFACE CROSSING
                CALL collect_char
     &            (ipid,iplot,elt(iplot),ifa,eta(iplot),fre(iplot),isp,
     &             nsp,xplot(iplot),yplot(iplot),
     &             zplot(iplot),fplot(iplot),
     &             dx(iplot),dy(iplot),dz(iplot),df(iplot),
     &             ifapar,nchdim,nchara)
              ELSE
!               A LOST-AGAIN TRACEBACK DETECTED
!               PROCESSOR NUMBER
                recvchar(iplot)%NEPID=ifapar(ifa,elt(iplot))
                recvchar(iplot)%INE=ifapar(ifa+3,elt(iplot))
                recvchar(iplot)%KNE=eta(iplot)
                recvchar(iplot)%IFR=fre(iplot)
              ENDIF
!             EXITING LOOP ON ISP
              EXIT
            ENDIF
!
!-----------------------------------------------------------------------
!           EXIT THROUGH TOP OR BOTTOM, OR LIQUID OR SOLID BOUNDARIES
!-----------------------------------------------------------------------
!
            dxp = dx(iplot)
            dyp = dy(iplot)
            i1  = ikle2(elt(iplot),ifa)
            i2  = ikle2(elt(iplot),isui(ifa))
            dx1 = x(i2) - x(i1)
            dy1 = y(i2) - y(i1)
!
            IF(iel.EQ.-1) THEN
!
!-----------------------------------------------------------------------
!  HERE: THE EXIT FACE IS A SOLID BOUNDARY
!  SETS SHP TO 0, END OF TRACING BACK (PROVISIONAL !!!!!!)
!-----------------------------------------------------------------------
!
              shp(1,iplot) = 0.d0
              shp(2,iplot) = 0.d0
              shp(3,iplot) = 0.d0
              elt(iplot) = - sens * elt(iplot)
              EXIT   ! LOOP ON ISP
!
            ENDIF
!
!-----------------------------------------------------------------------
!  HERE: THE EXIT FACE IS A LIQUID BOUNDARY
!  ENDS TRACING BACK (SIGN OF ELT)
!-----------------------------------------------------------------------
!
            denom = dxp*dy1-dyp*dx1
            IF(abs(denom).GT.1.d-12) THEN
              a1  = (dxp*(yp-y(i1))-dyp*(xp-x(i1))) / denom
            ELSE
              a1  = 0.d0
            ENDIF
            IF(a1.GT.1.d0) a1 = 1.d0
            IF(a1.LT.0.d0) a1 = 0.d0
            shp(ifa            ,iplot) = 1.d0 - a1
            shp(isui(ifa)      ,iplot) = a1
            shp(isui(isui(ifa)),iplot) = 0.d0
            xplot(iplot) = x(i1) + a1 * dx1
            yplot(iplot) = y(i1) + a1 * dy1
            IF(abs(dxp).GT.abs(dyp)) THEN
              a1 = (xp-xplot(iplot))/dxp
            ELSE
              a1 = (yp-yplot(iplot))/dyp
            ENDIF
            zplot(iplot) = zp - a1*dz(iplot)
            shz(iplot) = (zplot(iplot)-zstar(iet))
     &                       / (zstar(iet+1)-zstar(iet))
            fplot(iplot) = fp - a1*df(iplot)
            shf(iplot) = (fplot(iplot)-freq(ifr))
     &                       / (freq(ifr+1)-freq(ifr))
            elt(iplot) = - sens * elt(iplot)
!           END OF TRACING BACK
            EXIT   ! LOOP ON ISP
!
          ELSEIF(ifa.LE.5) THEN
!
!         IEL = IBOR(IEL,IFA,IET)
          iel = ibor(iel,ifa,1)
!
!-----------------------------------------------------------------------
!  HERE: THE EXIT FACE OF THE PRISM IS A TRIANGULAR FACE IN Z
!  =====================================================================
!-----------------------------------------------------------------------
!
            ifa = ifa - 4
!
            IF(iel.EQ.1) THEN
!
!-----------------------------------------------------------------------
!  HERE: THE EXIT FACE IS AN INTERIOR FACE
!  MOVES TO THE ADJACENT ELEMENT
!-----------------------------------------------------------------------
!
              eta(iplot) = iet + ifa + ifa - 1
              IF(eta(iplot).EQ.nplan+1) THEN
                eta(iplot)=1
                zp=zp-zstar(nplan+1)
                zplot(iplot)=zp
              ENDIF
              IF(eta(iplot).EQ.0) THEN
                eta(iplot) = nplan
                zp=zp+zstar(nplan+1)
                zplot(iplot)=zp
              ENDIF
              shz(iplot) = (zp-zstar(eta(iplot)))
     &                   / (zstar(eta(iplot)+1)-zstar(eta(iplot)))
!
              IF(add) THEN
                recvchar(iplot)%KNE=eta(iplot)
                recvchar(iplot)%ZP=zp
              ENDIF
!
              GO TO 50
!
            ELSE
              WRITE(lu,*) 'PROBLEM IN SCHAR41_PER_4D',iel,iplot
              CALL plante(1)
              stop
            ENDIF
!
          ELSE
!
!-----------------------------------------------------------------------
!  HERE: THE EXIT FACE OF THE PRISM IS A TRIANGULAR FACE FREQ
!  =====================================================================
!-----------------------------------------------------------------------
!
!           IBOR IS NOT REALLY BUILT FOR IFA = 6 OR 7 BUT IS ALWAYS 1.
!           IEL = IBOR(IEL,IFA,1)
            iel=1
            ifa = ifa - 6
!
            IF(ifa.EQ.1.AND.ifr.EQ.nf-1) iel=-1
            IF(ifa.EQ.0.AND.ifr.EQ.1)    iel=-1
            IF(iel.EQ.1) THEN
!
!-----------------------------------------------------------------------
!  HERE: THE EXIT FACE IS AN INTERIOR FACE
!  MOVES TO THE ADJACENT ELEMENT
!-----------------------------------------------------------------------
!
              fre(iplot) = ifr + ifa + ifa - 1
              shf(iplot) = (fp-freq(fre(iplot)))
     &                  / (freq(fre(iplot)+1)-freq(fre(iplot)))
!
              IF(add) THEN
                recvchar(iplot)%IFR=fre(iplot)
              ENDIF
!
              GOTO 50
!
            ELSE
!
!-----------------------------------------------------------------------
!  HERE: THE EXIT FACE IS THE MIN OR MAX FREQUENCY
!  PROJECTS THE REMAINING PART ON THE BOUNDARY AND CONTINUES
!-----------------------------------------------------------------------
!
              fplot(iplot)=freq(ifr+ifa)
              df(iplot)=0.d0
              shf(iplot)=ifa
              iso = isoh +isot
              IF(iso.NE.0) THEN
                GO TO 50
              ENDIF
!
            ENDIF
!
          ENDIF
!
!       IF(ISO.NE.0) THEN
        ENDIF
!
        ENDDO
      ENDDO
!
      IF(.NOT.add) THEN
        IF(ncsize.GT.1) nspmax=p_max(nspmax)
        WRITE(lu,*) 'NUMBER OF SUB-ITERATIONS :',nspmax
      ENDIF
!
!     RESTORING ORIGINAL ETA1
!
      eta1(nplan)=nplan+1
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE schar41_per_4d
!
!-----------------------------------------------------------------------
! CHAR11 MODIFIED FOR INITIAL COLLECTING OF THE LOST CHARACTERISTICS
! I.E. THE ONES CROSSING INTERFACE PARTITIONS IN THE PARALLEL CASE
! IFAPAR  :: DELIVERS LOCAL ELEMENT NUMBER AND THE PARTITION NR THERE
!            WHEN CROSSING THE INTERFACE VIA A HALO ELEMENT FACE
!-----------------------------------------------------------------------
! JAJ MODIFIED WED JUL 16 18:24:08 CEST 2008
!
!                       ******************
                        SUBROUTINE schar11
!                       ******************
!
     &(u,v,dt,nrk,x,y,ikle,ifabor,xplot,yplot,dx,dy,shp,elt,
     & nplot,nelem,nelmax,surdet,sens,
     & ifapar,nchdim,nchara,add)
!
!***********************************************************************
! BIEF VERSION 6.2           24/04/97    J-M JANIN (LNH) 30 87 72 84
!
!***********************************************************************
!
!  FONCTION :
!
!     REMONTEE OU DESCENTE
!     DES COURBES CARACTERISTIQUES
!     SUR DES QUADRILATERES P1
!     DANS L'INTERVALLE DE TEMPS DT
!     AVEC UNE DISCRETISATION ELEMENTS FINIS
!
!
!  DISCRETISATION :
!
!     LE DOMAINE EST APPROCHE PAR UNE DISCRETISATION ELEMENTS FINIS
!     UNE APPROXIMATION LOCALE EST DEFINIE POUR LE VECTEUR VITESSE :
!     LA VALEUR EN UN POINT D'UN ELEMENT NE DEPEND QUE DES VALEURS
!     AUX NOEUDS DE CET ELEMENT
!
!
!  RESTRICTIONS ET HYPOTHESES :
!
!     LE CHAMP CONVECTEUR U EST SUPPOSE INDEPENDANT DU TEMPS
!     LE DERIVANT EST SUPPOSE PONCTUEL DONC NON DISPERSIF
!
!-----------------------------------------------------------------------
!                             ARGUMENTS
! .________________.____.______________________________________________.
! |      NOM       |MODE|                   ROLE                       |
! |________________|____|______________________________________________|
! |    U,V         | -->| COMPOSANTE DE LA VITESSE DU CONVECTEUR       |
! |    DT          | -->| PAS DE TEMPS.
! |    NRK         | -->| NOMBRE DE SOUS-PAS DE RUNGE-KUTTA.
! |    X,Y         | -->| COORDONNEES DES POINTS DU MAILLAGE.
! |    IKLE        | -->| TRANSITION ENTRE LES NUMEROTATIONS LOCALE
! |                |    | ET GLOBALE.
! |    IFABOR      | -->| NUMEROS DES ELEMENTS AYANT UNE FACE COMMUNE
! |                |    | AVEC L'ELEMENT .  SI IFABOR<0 OU NUL
! |                |    | ON A UNE FACE LIQUIDE,SOLIDE,OU PERIODIQUE
! |  XPLOT,YPLOT   |<-->| POSITIONS SUCCESSIVES DES DERIVANTS.
! |    DX,DY       | -- | STOCKAGE DES SOUS-PAS . |
! |    SHP         |<-->| COORDONNEES BARYCENTRIQUES 2D AU PIED DES
! |                |    | COURBES CARACTERISTIQUES.
! |    ELT         |<-->| NUMEROS DES ELEMENTS 2D AU PIED DES COURBES
! |                |    | CARACTERISTIQUES.
! |    NPLOT       | -->| NOMBRE DE DERIVANTS.
! |    NELEM       | -->| NOMBRE D'ELEMENTS.
! |    NELMAX      | -->| NOMBRE MAXIMAL D'ELEMENTS DANS LE MAILLAGE 2D
! |    SURDET      | -->| VARIABLE UTILISEE PAR LA TRANSFORMEE ISOPARAM.
! |    SENS        | -->| DESCENTE OU REMONTEE DES CARACTERISTIQUES.
! |________________|____|______________________________________________|
!  MODE: -->(DONNEE NON MODIFIEE),<--(RESULTAT),<-->(DONNEE MODIFIEE)
!-----------------------------------------------------------------------
!     - APPELE PAR : CARACT , DERIVE , DERLAG
!     - PROGRAMMES APPELES : NEANT
!
!***********************************************************************
!
      USE bief
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)    :: SENS,NCHDIM
      INTEGER         , INTENT(IN)    :: NELEM,NELMAX,NPLOT,NRK
      INTEGER         , INTENT(IN)    :: IKLE(nelmax,3),IFABOR(nelmax,3)
      INTEGER         , INTENT(INOUT) :: ELT(*),NCHARA
      DOUBLE PRECISION, INTENT(IN)    :: U(*),V(*),SURDET(nelem)
      DOUBLE PRECISION, INTENT(INOUT) :: XPLOT(*),YPLOT(*)
      DOUBLE PRECISION, INTENT(INOUT) :: SHP(3,*)
      DOUBLE PRECISION, INTENT(IN)    :: DT
      DOUBLE PRECISION, INTENT(IN)    :: X(*),Y(*)
      DOUBLE PRECISION, INTENT(INOUT) :: DX(*),DY(*)
      INTEGER, INTENT(IN)             :: IFAPAR(6,*)
      LOGICAL, INTENT(IN)             :: ADD
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER :: IPLOT,ISP,I1,I2,I3,IEL,ISO,IFA,ISUI(3),ISUI2(3),ISPDONE
      INTEGER IPROC,ILOC,NSP
      DOUBLE PRECISION PAS,A1,DX1,DY1,DXP,DYP,XP,YP,DENOM
!
      INTRINSIC int,max,min,sqrt
!
      parameter( isui = (/ 2 , 3 , 1 /) )
      parameter( isui2 = (/ 3 , 1 , 2 /) )
      DOUBLE PRECISION, PARAMETER :: EPSILO = -1.d-6
!
!-----------------------------------------------------------------------
!     FOR EVERY POINT
!-----------------------------------------------------------------------
!
      DO iplot=1,nplot
!
        IF(add) THEN
!
          xplot(iplot)   = recvchar(iplot)%XP
          yplot(iplot)   = recvchar(iplot)%YP
          dx(iplot)      = recvchar(iplot)%DX
          dy(iplot)      = recvchar(iplot)%DY
          elt(iplot)     = recvchar(iplot)%INE
          nsp            = recvchar(iplot)%NSP ! R-K STEPS TO BE FULLFILLED
          ispdone        = recvchar(iplot)%ISP ! R-K STEPS ALREADY DONE
          iel = elt(iplot)
          xp  = xplot(iplot)
          yp  = yplot(iplot)
          i1 = ikle(iel,1)
          i2 = ikle(iel,2)
          i3 = ikle(iel,3)
          shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &           -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
          shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &           -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
          shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &           -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
!         ASSUME TO BE LOCALISED IT WILL BE SET OTHERWISE IF LOST-AGAIN
          recvchar(iplot)%NEPID=-1
!
        ELSE
          iel = elt(iplot)
!         POINTS WITH IEL=0 ARE TREATED SO THAT THE FINAL INTERPOLATION
!         GIVES 0., AND WE SKIP TO NEXT POINT IPLOT (CYCLE)
!         THIS WILL NOT INTERFERE WITH ELT(IPLOT)=0 GIVEN ON LIQUID BOUNDARIES
!         BY ARRAY IFABOR, THAT MAY HAPPEN LATER
          IF(iel.EQ.0) THEN
            elt(iplot)=1
            shp(1,iplot)=0.d0
            shp(2,iplot)=0.d0
            shp(3,iplot)=0.d0
            cycle
          ENDIF
          i1 = ikle(iel,1)
          i2 = ikle(iel,2)
          i3 = ikle(iel,3)
          dxp = u(i1)*shp(1,iplot)+u(i2)*shp(2,iplot)
     &                            +u(i3)*shp(3,iplot)
          dyp = v(i1)*shp(1,iplot)+v(i2)*shp(2,iplot)
     &                            +v(i3)*shp(3,iplot)
          nsp=max(1,int(nrk*dt*sqrt((dxp**2+dyp**2)*surdet(iel))))
          ispdone=1
        ENDIF
!
        pas = sens * dt / nsp
!
!       LOOP ON RUNGE-KUTTA SUB-STEPS
!
!       COMPILER MUST DO NOTHING IF ISPDONE>NSP
!       IN MODE "ADD", ISP = ISPDONE HAS NOT BEEN FULLY DONE
!       IT IS RESTARTED HERE
!
        DO isp=ispdone,nsp
!
!-----------------------------------------------------------------------
!       LOCALISING THE ARRIVAL POINT
!-----------------------------------------------------------------------
!
!       IN MODE "ADD" ITERATIONS ALREADY DONE ARE SKIPPED AND
!                     CHARACTERISTICS GONE IN ANOTHER SUB-DOMAIN SKIPPED
!
        IF(add) THEN
          IF(isp.EQ.ispdone) GO TO 50
          IF(recvchar(iplot)%NEPID.NE.-1) cycle
        ENDIF
!
        iel = elt(iplot)
        i1 = ikle(iel,1)
        i2 = ikle(iel,2)
        i3 = ikle(iel,3)
!
        dx(iplot) = ( u(i1)*shp(1,iplot)
     &              + u(i2)*shp(2,iplot)
     &              + u(i3)*shp(3,iplot) ) * pas
        dy(iplot) = ( v(i1)*shp(1,iplot)
     &              + v(i2)*shp(2,iplot)
     &              + v(i3)*shp(3,iplot) ) * pas
        xp = xplot(iplot) + dx(iplot)
        yp = yplot(iplot) + dy(iplot)
!
        shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                 -(y(i3)-y(i2))*(xp-x(i2))) * surdet(iel)
        shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                 -(y(i1)-y(i3))*(xp-x(i3))) * surdet(iel)
        shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                 -(y(i2)-y(i1))*(xp-x(i1))) * surdet(iel)
!
        xplot(iplot) = xp
        yplot(iplot) = yp
!
        IF(add) THEN
!         CONTINUOUS SETTING OF THE REACHED POSITION FOR IPLOT
!         AND THE NUMBER OF STEPS DONE ALREADY
          recvchar(iplot)%XP=xplot(iplot)
          recvchar(iplot)%YP=yplot(iplot)
          recvchar(iplot)%DX=dx(iplot)
          recvchar(iplot)%DY=dy(iplot)
          recvchar(iplot)%INE=elt(iplot)
          recvchar(iplot)%ISP=isp
        ENDIF
!
!-----------------------------------------------------------------------
!       TEST: IS THE PATHLINE WENT OUT THE ORIGINAL ELEMENT
!-----------------------------------------------------------------------
!
50      CONTINUE
!
        iso = 0
        IF(shp(1,iplot).LT.epsilo) iso = 1
        IF(shp(2,iplot).LT.epsilo) iso = iso + 2
        IF(shp(3,iplot).LT.epsilo) iso = iso + 4
!
        IF(iso.NE.0) THEN
!
!-----------------------------------------------------------------------
!         HERE WE ARE OUT OF THE ELEMENT
!-----------------------------------------------------------------------
!
          iel = elt(iplot)
          xp = xplot(iplot)
          yp = yplot(iplot)
!
!         THE 3 LINES FORMING THE TRIANGLE CUT THE PLANE INTO 7
!         ZONES, NUMBERED FROM 0 (INSIDE THE TRIANGLE) TO 6
!         ISO IS THE NUMBER. FOR ISO =1,2,4, THERE IS NO AMBIGUITY
!         AS TO THE EDGE CROSSED. FOR ISO = 3, IT CAN BE EDGE 2
!         OR 3, FOR ISO = 5 IT CAN BE EDGE 1 OR 2, FOR ISO = 6 IT
!         CAN BE EDGE 1 OR 3.
!         FOR CASES 3, 5 AND 6, AN INNER PRODUCT SHOWS IF THE DIRECTION
!         OF THE DISPLACEMENT (DX,DY) IS ON THE RIGHT OR ON THE LEFT
!         OF THE INTERSECTION BETWEEN THE TWO EDGES, SO IT GIVES
!         THE REAL EDGE THAT HAS BEEN CROSSED
!
          IF(iso.EQ.1) THEN
            ifa = 2
          ELSEIF (iso.EQ.2) THEN
            ifa = 3
          ELSEIF (iso.EQ.4) THEN
            ifa = 1
          ELSEIF (iso.EQ.3) THEN
            ifa = 2
            IF(dx(iplot)*(y(ikle(iel,3))-yp).LT.
     &         dy(iplot)*(x(ikle(iel,3))-xp)) ifa = 3
          ELSEIF (iso.EQ.6) THEN
            ifa = 3
            IF (dx(iplot)*(y(ikle(iel,1))-yp).LT.
     &          dy(iplot)*(x(ikle(iel,1))-xp)) ifa = 1
          ELSE
!           HERE CASE ISO=5
            ifa = 1
            IF(dx(iplot)*(y(ikle(iel,2))-yp).LT.
     &         dy(iplot)*(x(ikle(iel,2))-xp)) ifa = 2
          ENDIF
!
          iel = ifabor(iel,ifa)
!
          IF(iel.GT.0) THEN
!
!-----------------------------------------------------------------------
!           HERE WE ARRIVE IN ANOTHER ELEMENT
!-----------------------------------------------------------------------
!
            i1 = ikle(iel,1)
            i2 = ikle(iel,2)
            i3 = ikle(iel,3)
!
            elt(iplot) = iel
            shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                     -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
            shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                     -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
            shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                     -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
!
            GOTO 50
!
          ENDIF
!
!-----------------------------------------------------------------------
!         HERE WE PASS TO NEIGHBOUR SUBDOMAIN AND COLLECT DATA
!-----------------------------------------------------------------------
!
          IF(iel.EQ.-2) THEN
            IF(add) THEN
!             A LOST-AGAIN TRACEBACK DETECTED, ALREADY HERE
!             SET THE IMPLANTING PARAMETERS
              iproc=ifapar(ifa  ,elt(iplot))
              iloc =ifapar(ifa+3,elt(iplot))
!             ANOTHER ONE AS IPID, MEANS ALSO NOT LOCALISED
              recvchar(iplot)%NEPID=iproc
              recvchar(iplot)%INE=iloc
            ELSE
              CALL collect_char(ipid,iplot,elt(iplot),ifa,0,0,isp,nsp,
     &                          xplot(iplot),yplot(iplot),0.d0,0.d0,
     &                          dx(iplot),dy(iplot),0.d0,0.d0,
     &                          ifapar,nchdim,nchara)
            ENDIF
!           EXITING LOOP ON ISP
            EXIT
          ENDIF
!
!-----------------------------------------------------------------------
!         SPECIAL TREATMENT FOR SOLID OR LIQUID BOUNDARIES
!-----------------------------------------------------------------------
!
          dxp = dx(iplot)
          dyp = dy(iplot)
          i1  = ikle(elt(iplot),ifa)
          i2  = ikle(elt(iplot),isui(ifa))
          dx1 = x(i2) - x(i1)
          dy1 = y(i2) - y(i1)
!
          IF(iel.EQ.-1) THEN
!
!-----------------------------------------------------------------------
!           HERE SOLID BOUNDARY, VELOCITY IS PROJECTED ON THE BOUNDARY
!           AND WE GO ON
!-----------------------------------------------------------------------
!
!           HERE A1 IS THE PARAMETRIC COORDINATE OF THE PROJECTED
!           DISPLACEMENT ON SEGMENT I1----I2
!
            a1 = (dxp*dx1 + dyp*dy1) / (dx1**2 + dy1**2)
!
!           THE TOTAL DISPLACEMENT IS PROJECTED HERE, NOT THE REMAINING
!           PART, BUT ONLY THE DIRECTION WILL BE USED
            dx(iplot) = a1 * dx1
            dy(iplot) = a1 * dy1
!
!           NOW A1 IS THE PARAMETRIC COORDINATE ON SEGMENT I1----I2
!           OF THE POSITION OF THE ARRIVAL POINT, I.E. INTERSECTION
!           + REMAINING DISPLACEMENT PROJECTED ON THE SEGMENT
!           ITS VALUE MAY BE OUTSIDE THE RANGE (0,1). THE VALUE OF A1
!           SIMPLIFIES INTO THE FOLLOWING FORMULA, BECAUSE IT IS
!           SIMPLY VECTOR I1----P PROJECTED ON SEGMENT I1----I2
!
            a1 = ((xp-x(i1))*dx1+(yp-y(i1))*dy1)/(dx1**2+dy1**2)
            shp(      ifa ,iplot) = 1.d0 - a1
            shp( isui(ifa),iplot) = a1
            shp(isui2(ifa),iplot) = 0.d0
            xplot(iplot) = x(i1) + a1 * dx1
            yplot(iplot) = y(i1) + a1 * dy1
            IF(add) THEN
              recvchar(iplot)%XP=xplot(iplot)
              recvchar(iplot)%YP=yplot(iplot)
              recvchar(iplot)%DX=dx(iplot)
              recvchar(iplot)%DY=dy(iplot)
            ENDIF
!
            GOTO 50
!
          ELSEIF(iel.EQ.0) THEN
!
!------------------------------------------------------------------------
!           HERE WE HAVE A LIQUID BOUNDARY, THE CHARACTERISTIC IS STOPPED
!------------------------------------------------------------------------
!
            denom = dxp*dy1-dyp*dx1
            IF(abs(denom).GT.1.d-12) THEN
              a1  = (dxp*(yp-y(i1))-dyp*(xp-x(i1))) / denom
            ELSE
              a1  = 0.d0
            ENDIF
            a1 = max(min(a1,1.d0),0.d0)
            shp(      ifa ,iplot) = 1.d0 - a1
            shp( isui(ifa),iplot) = a1
            shp(isui2(ifa),iplot) = 0.d0
            xplot(iplot) = x(i1) + a1 * dx1
            yplot(iplot) = y(i1) + a1 * dy1
!           THIS IS A MARKER FOR PARTICLES EXITING A DOMAIN
!           SENS=-1 FOR BACKWARD CHARACTERISTICS
            elt(iplot) = - sens * elt(iplot)
            EXIT
!
          ELSE
!
            WRITE(lu,*) 'UNEXPECTED CASE IN SCHAR11'
            WRITE(lu,*) 'IEL=',iel
            CALL plante(1)
            stop
!
          ENDIF
!
        ENDIF
!
        ENDDO
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE schar11
!                       **********************
                        SUBROUTINE schar11_sto
!                       **********************
!
     &(u,v,dt,nrk,x,y,ikle,ifabor,xplot,yplot,dx,dy,shp,elt,
     & nplot,npoin,nelem,nelmax,surdet,sens,
     & ifapar,nchdim,nchara,add,ielm,visc,stocha)
!
!***********************************************************************
! BIEF VERSION 6.2           24/04/97    J-M JANIN (LNH) 30 87 72 84
!
!***********************************************************************
!
!  FONCTION :
!
!     REMONTEE OU DESCENTE
!     DES COURBES CARACTERISTIQUES
!     SUR DES QUADRILATERES P1
!     DANS L'INTERVALLE DE TEMPS DT
!     AVEC UNE DISCRETISATION ELEMENTS FINIS
!
!
!  DISCRETISATION :
!
!     LE DOMAINE EST APPROCHE PAR UNE DISCRETISATION ELEMENTS FINIS
!     UNE APPROXIMATION LOCALE EST DEFINIE POUR LE VECTEUR VITESSE :
!     LA VALEUR EN UN POINT D'UN ELEMENT NE DEPEND QUE DES VALEURS
!     AUX NOEUDS DE CET ELEMENT
!
!
!  RESTRICTIONS ET HYPOTHESES :
!
!     LE CHAMP CONVECTEUR U EST SUPPOSE INDEPENDANT DU TEMPS
!     LE DERIVANT EST SUPPOSE PONCTUEL DONC NON DISPERSIF
!
!-----------------------------------------------------------------------
!                             ARGUMENTS
! .________________.____.______________________________________________.
! |      NOM       |MODE|                   ROLE                       |
! |________________|____|______________________________________________|
! |    U,V         | -->| COMPOSANTE DE LA VITESSE DU CONVECTEUR       |
! |    DT          | -->| PAS DE TEMPS.
! |    NRK         | -->| NOMBRE DE SOUS-PAS DE RUNGE-KUTTA.
! |    X,Y         | -->| COORDONNEES DES POINTS DU MAILLAGE.
! |    IKLE        | -->| TRANSITION ENTRE LES NUMEROTATIONS LOCALE
! |                |    | ET GLOBALE.
! |    IFABOR      | -->| NUMEROS DES ELEMENTS AYANT UNE FACE COMMUNE
! |                |    | AVEC L'ELEMENT .  SI IFABOR<0 OU NUL
! |                |    | ON A UNE FACE LIQUIDE,SOLIDE,OU PERIODIQUE
! |  XPLOT,YPLOT   |<-->| POSITIONS SUCCESSIVES DES DERIVANTS.
! |    DX,DY       | -- | STOCKAGE DES SOUS-PAS . |
! |    SHP         |<-->| COORDONNEES BARYCENTRIQUES 2D AU PIED DES
! |                |    | COURBES CARACTERISTIQUES.
! |    ELT         |<-->| NUMEROS DES ELEMENTS 2D AU PIED DES COURBES
! |                |    | CARACTERISTIQUES.
! |    NPLOT       | -->| NOMBRE DE DERIVANTS.
! |    NPOIN       | -->| NOMBRE DE POINTS DU MAILLAGE.
! |    NELEM       | -->| NOMBRE D'ELEMENTS.
! |    NELMAX      | -->| NOMBRE MAXIMAL D'ELEMENTS DANS LE MAILLAGE 2D
! |    SENS        | -->| -1: BACKWARD CHARACTERISTICS 1: FORWARD
! |    STOCHA      | -->| STOCHASTIC DIFFUSION MODEL
! |    SURDET      | -->| VARIABLE UTILISEE PAR LA TRANSFORMEE ISOPARAM.
! |________________|____|______________________________________________|
!  MODE: -->(DONNEE NON MODIFIEE),<--(RESULTAT),<-->(DONNEE MODIFIEE)
!-----------------------------------------------------------------------
!     - APPELE PAR : CARACT , DERIVE , DERLAG
!     - PROGRAMMES APPELES : NEANT
!
!***********************************************************************
!
      USE bief
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)    :: SENS,NCHDIM,IELM,STOCHA
      INTEGER         , INTENT(IN)    :: NPOIN,NELEM,NELMAX,NPLOT,NRK
      INTEGER         , INTENT(IN)    :: IKLE(nelmax,*),IFABOR(nelmax,3)
      INTEGER         , INTENT(INOUT) :: ELT(*),NCHARA
      DOUBLE PRECISION, INTENT(IN)    :: U(npoin),V(npoin),SURDET(nelem)
      DOUBLE PRECISION, INTENT(INOUT) :: XPLOT(*),YPLOT(*)
      DOUBLE PRECISION, INTENT(INOUT) :: SHP(3,*)
      DOUBLE PRECISION, INTENT(IN)    :: DT
      DOUBLE PRECISION, INTENT(IN)    :: X(npoin),Y(npoin),VISC(npoin)
      DOUBLE PRECISION, INTENT(INOUT) :: DX(*),DY(*)
      INTEGER, INTENT(IN)             :: IFAPAR(6,*)
      LOGICAL, INTENT(IN)             :: ADD
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IPLOT,ISP,I1,I2,I3,I4,I5,I6
      INTEGER :: IEL,ISO,IFA,ISUI(3),ISUI2(3),ISPDONE
      INTEGER IPROC,ILOC,NSP
      DOUBLE PRECISION PAS,A1,DX1,DY1,DXP,DYP,XP,YP,DENOM
      DOUBLE PRECISION SHP11,SHP12,SHP14
      DOUBLE PRECISION SHP22,SHP23,SHP24
      DOUBLE PRECISION SHP33,SHP31,SHP34
!     FOR STOCHASTIC DIFFUSION
      DOUBLE PRECISION RAND1,RAND2,A,C,D,E,DIFF_X,DIFF_Y,DEUXPI
!
      INTRINSIC int,max,min,sqrt,acos
!
      parameter( isui = (/ 2 , 3 , 1 /) )
      parameter( isui2 = (/ 3 , 1 , 2 /) )
      DOUBLE PRECISION, PARAMETER :: EPSILO = -1.d-6
!
!-----------------------------------------------------------------------
!     FOR EVERY POINT
!-----------------------------------------------------------------------
!
      deuxpi=2.d0*acos(-1.d0)
!
      DO iplot=1,nplot
!
        IF(add) THEN
!
          xplot(iplot)   = recvchar(iplot)%XP
          yplot(iplot)   = recvchar(iplot)%YP
          dx(iplot)      = recvchar(iplot)%DX
          dy(iplot)      = recvchar(iplot)%DY
          elt(iplot)     = recvchar(iplot)%INE
          nsp            = recvchar(iplot)%NSP ! R-K STEPS TO BE FULLFILLED
          ispdone        = recvchar(iplot)%ISP ! R-K STEPS ALREADY DONE
          iel = elt(iplot)
          xp  = xplot(iplot)
          yp  = yplot(iplot)
          i1 = ikle(iel,1)
          i2 = ikle(iel,2)
          i3 = ikle(iel,3)
          shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &           -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
          shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &           -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
          shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &           -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
!         ASSUME TO BE LOCALISED IT WILL BE SET OTHERWISE IF LOST-AGAIN
          recvchar(iplot)%NEPID=-1
!
        ELSE
          iel = elt(iplot)
!         POINTS WITH IEL=0 ARE TREATED SO THAT THE FINAL INTERPOLATION
!         GIVES 0., AND WE SKIP TO NEXT POINT IPLOT (CYCLE)
!         THIS WILL NOT INTERFERE WITH ELT(IPLOT)=0 GIVEN ON LIQUID BOUNDARIES
!         BY ARRAY IFABOR, THAT MAY HAPPEN LATER
          IF(iel.EQ.0) THEN
            elt(iplot)=1
            shp(1,iplot)=0.d0
            shp(2,iplot)=0.d0
            shp(3,iplot)=0.d0
            cycle
          ENDIF
          i1 = ikle(iel,1)
          i2 = ikle(iel,2)
          i3 = ikle(iel,3)
!         HERE WITHOUT CONSIDERING STOCHASTIC DIFFUSION
!         NOR QUASI-BUBBLE OR QUADRATIC DISCRETISATION
!         THIS IS JUST FOR COMPUTING NSP
          dxp = u(i1)*shp(1,iplot)+u(i2)*shp(2,iplot)
     &                            +u(i3)*shp(3,iplot)
          dyp = v(i1)*shp(1,iplot)+v(i2)*shp(2,iplot)
     &                            +v(i3)*shp(3,iplot)
          nsp=max(1,int(nrk*dt*sqrt((dxp**2+dyp**2)*surdet(iel))))
          ispdone=1
        ENDIF
!
        pas = sens * dt / nsp
!
!       LOOP ON RUNGE-KUTTA SUB-STEPS
!
!       COMPILER MUST DO NOTHING IF ISPDONE>NSP
!       IN MODE "ADD", ISP = ISPDONE HAS NOT BEEN FULLY DONE
!       IT IS RESTARTED HERE
!
        DO isp=ispdone,nsp
!
!-----------------------------------------------------------------------
!       LOCALISING THE ARRIVAL POINT
!-----------------------------------------------------------------------
!
!       IN MODE "ADD" ITERATIONS ALREADY DONE ARE SKIPPED AND
!                     CHARACTERISTICS GONE IN ANOTHER SUB-DOMAIN SKIPPED
!
        IF(add) THEN
          IF(isp.EQ.ispdone) GO TO 50
          IF(recvchar(iplot)%NEPID.NE.-1) cycle
        ENDIF
!
        iel = elt(iplot)
        i1 = ikle(iel,1)
        i2 = ikle(iel,2)
        i3 = ikle(iel,3)
!
        IF(ielm.EQ.11) THEN
!
          dx(iplot) = ( u(i1)*shp(1,iplot)
     &                + u(i2)*shp(2,iplot)
     &                + u(i3)*shp(3,iplot) ) * pas
          dy(iplot) = ( v(i1)*shp(1,iplot)
     &                + v(i2)*shp(2,iplot)
     &                + v(i3)*shp(3,iplot) ) * pas
!
        ELSEIF(ielm.EQ.12) THEN
!
          i4 = ikle(iel,4)
          shp11=shp(1,iplot)-shp(3,iplot)
          shp12=shp(2,iplot)-shp(3,iplot)
          shp14=3.d0*shp(3,iplot)
          shp22=shp(2,iplot)-shp(1,iplot)
          shp23=shp(3,iplot)-shp(1,iplot)
          shp24=3.d0*shp(1,iplot)
          shp33=shp(3,iplot)-shp(2,iplot)
          shp31=shp(1,iplot)-shp(2,iplot)
          shp34=3.d0*shp(2,iplot)
          IF(     shp11.GT.     2.d0*epsilo .AND.
     &            shp11.LT.1.d0-4.d0*epsilo .AND.
     &            shp12.GT.     2.d0*epsilo .AND.
     &            shp12.LT.1.d0-4.d0*epsilo .AND.
     &            shp14.LT.1.d0-4.d0*epsilo ) THEN
          dx(iplot) = ( u(i1) * shp11
     &                + u(i2) * shp12
     &                + u(i4) * shp14 ) * pas
          dy(iplot) = ( v(i1) * shp11
     &                + v(i2) * shp12
     &                + v(i4) * shp14 ) * pas
          ELSEIF( shp22.GT.     2.d0*epsilo .AND.
     &            shp22.LT.1.d0-4.d0*epsilo .AND.
     &            shp23.GT.     2.d0*epsilo .AND.
     &            shp23.LT.1.d0-4.d0*epsilo .AND.
     &            shp24.LT.1.d0-4.d0*epsilo ) THEN
            dx(iplot) = ( u(i2) * shp22
     &                  + u(i3) * shp23
     &                  + u(i4) * shp24 ) * pas
            dy(iplot) = ( v(i2) * shp22
     &                  + v(i3) * shp23
     &                  + v(i4) * shp24 ) * pas
          ELSEIF( shp33.GT.     2.d0*epsilo .AND.
     &            shp33.LT.1.d0-4.d0*epsilo .AND.
     &            shp31.GT.     2.d0*epsilo .AND.
     &            shp31.LT.1.d0-4.d0*epsilo .AND.
     &            shp34.LT.1.d0-4.d0*epsilo ) THEN
            dx(iplot) = ( u(i3) * shp33
     &                  + u(i1) * shp31
     &                  + u(i4) * shp34 ) * pas
            dy(iplot) = ( v(i3) * shp33
     &                  + v(i1) * shp31
     &                  + v(i4) * shp34 ) * pas
          ELSE
            WRITE(lu,*) 'SCHAR11_STO: POINT ',iplot
            WRITE(lu,*) 'NOT IN ELEMENT ',elt(iplot)
            WRITE(lu,*) 'SHP(1,IPLOT)=',shp(1,iplot)
            WRITE(lu,*) 'SHP(2,IPLOT)=',shp(2,iplot)
            WRITE(lu,*) 'SHP(3,IPLOT)=',shp(3,iplot)
            WRITE(lu,*) 'EPSILO=',epsilo,' IPID=',ipid
            CALL plante(1)
            stop
          ENDIF
!
        ELSEIF(ielm.EQ.13) THEN
!
          i4 = ikle(iel,4)
          i5 = ikle(iel,5)
          i6 = ikle(iel,6)
          dx(iplot) = ( u(i1)*(2.d0*shp(1,iplot)-1.d0)*shp(1,iplot)
     &                + u(i2)*(2.d0*shp(2,iplot)-1.d0)*shp(2,iplot)
     &                + u(i3)*(2.d0*shp(3,iplot)-1.d0)*shp(3,iplot)
     &                + u(i4)*4.d0*shp(1,iplot)*shp(2,iplot)
     &                + u(i5)*4.d0*shp(2,iplot)*shp(3,iplot)
     &                + u(i6)*4.d0*shp(3,iplot)*shp(1,iplot)) * pas
          dy(iplot) = ( v(i1)*(2.d0*shp(1,iplot)-1.d0)*shp(1,iplot)
     &                + v(i2)*(2.d0*shp(2,iplot)-1.d0)*shp(2,iplot)
     &                + v(i3)*(2.d0*shp(3,iplot)-1.d0)*shp(3,iplot)
     &                + v(i4)*4.d0*shp(1,iplot)*shp(2,iplot)
     &                + v(i5)*4.d0*shp(2,iplot)*shp(3,iplot)
     &                + v(i6)*4.d0*shp(3,iplot)*shp(1,iplot)) * pas
!
        ELSE
!
          WRITE(lu,*) 'UNEXPECTED CASE IN SCHAR11_STO'
          WRITE(lu,*) 'IELM=',ielm
          CALL plante(1)
          stop
!
        ENDIF
!
!       STOCHASTIC DIFFUSION
!
        IF(stocha.EQ.1) THEN
!         COMPUTING LOCAL VISCOSITY
          a=max(visc(i1)*shp(1,iplot)
     &     +visc(i2)*shp(2,iplot)
     &     +visc(i3)*shp(3,iplot),0.d0)
!         DISPLACEMENT DUE TO RANDOM DIFFUSION
          CALL random_number(rand1)
          CALL random_number(rand2)
          c=sqrt(-2.d0*log(rand1))
          d=c*cos(deuxpi*rand2)
          e=c*sin(deuxpi*rand2)
          diff_x=d*sqrt(2.d0*a/0.72d0)
          diff_y=e*sqrt(2.d0*a/0.72d0)
          dx(iplot) = dx(iplot) + diff_x*sqrt(abs(pas))
          dy(iplot) = dy(iplot) + diff_y*sqrt(abs(pas))
        ENDIF
!
        xp = xplot(iplot) + dx(iplot)
        yp = yplot(iplot) + dy(iplot)
!
        shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                 -(y(i3)-y(i2))*(xp-x(i2))) * surdet(iel)
        shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                 -(y(i1)-y(i3))*(xp-x(i3))) * surdet(iel)
        shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                 -(y(i2)-y(i1))*(xp-x(i1))) * surdet(iel)
!
        xplot(iplot) = xp
        yplot(iplot) = yp
!
        IF(add) THEN
!         CONTINUOUS SETTING OF THE REACHED POSITION FOR IPLOT
!         AND THE NUMBER OF STEPS DONE ALREADY
          recvchar(iplot)%XP=xplot(iplot)
          recvchar(iplot)%YP=yplot(iplot)
          recvchar(iplot)%DX=dx(iplot)
          recvchar(iplot)%DY=dy(iplot)
          recvchar(iplot)%INE=elt(iplot)
          recvchar(iplot)%ISP=isp
        ENDIF
!
!-----------------------------------------------------------------------
!       TEST: IS THE PATHLINE WENT OUT THE ORIGINAL ELEMENT
!-----------------------------------------------------------------------
!
50      CONTINUE
!
        iso = 0
        IF(shp(1,iplot).LT.epsilo) iso = 1
        IF(shp(2,iplot).LT.epsilo) iso = iso + 2
        IF(shp(3,iplot).LT.epsilo) iso = iso + 4
!
        IF(iso.NE.0) THEN
!
!-----------------------------------------------------------------------
!         HERE WE ARE OUT OF THE ELEMENT
!-----------------------------------------------------------------------
!
          iel = elt(iplot)
          xp = xplot(iplot)
          yp = yplot(iplot)
!
!         THE 3 LINES FORMING THE TRIANGLE CUT THE PLANE INTO 7
!         ZONES, NUMBERED FROM 0 (INSIDE THE TRIANGLE) TO 6
!         ISO IS THE NUMBER. FOR ISO =1,2,4, THERE IS NO AMBIGUITY
!         AS TO THE EDGE CROSSED. FOR ISO = 3, IT CAN BE EDGE 2
!         OR 3, FOR ISO = 5 IT CAN BE EDGE 1 OR 2, FOR ISO = 6 IT
!         CAN BE EDGE 1 OR 3.
!         FOR CASES 3, 5 AND 6, AN INNER PRODUCT SHOWS IF THE DIRECTION
!         OF THE DISPLACEMENT (DX,DY) IS ON THE RIGHT OR ON THE LEFT
!         OF THE INTERSECTION BETWEEN THE TWO EDGES, SO IT GIVES
!         THE REAL EDGE THAT HAS BEEN CROSSED
!
          IF(iso.EQ.1) THEN
            ifa = 2
          ELSEIF (iso.EQ.2) THEN
            ifa = 3
          ELSEIF (iso.EQ.4) THEN
            ifa = 1
          ELSEIF (iso.EQ.3) THEN
            ifa = 2
            IF(dx(iplot)*(y(ikle(iel,3))-yp).LT.
     &         dy(iplot)*(x(ikle(iel,3))-xp)) ifa = 3
          ELSEIF (iso.EQ.6) THEN
            ifa = 3
            IF (dx(iplot)*(y(ikle(iel,1))-yp).LT.
     &          dy(iplot)*(x(ikle(iel,1))-xp)) ifa = 1
          ELSE
!           HERE CASE ISO=5
            ifa = 1
            IF(dx(iplot)*(y(ikle(iel,2))-yp).LT.
     &         dy(iplot)*(x(ikle(iel,2))-xp)) ifa = 2
          ENDIF
!
          iel = ifabor(iel,ifa)
!
          IF(iel.GT.0) THEN
!
!-----------------------------------------------------------------------
!           HERE WE ARRIVE IN ANOTHER ELEMENT
!-----------------------------------------------------------------------
!
            i1 = ikle(iel,1)
            i2 = ikle(iel,2)
            i3 = ikle(iel,3)
!
            elt(iplot) = iel
            shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                     -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
            shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                     -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
            shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                     -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
!
            GOTO 50
!
          ENDIF
!
!-----------------------------------------------------------------------
!         HERE WE PASS TO NEIGHBOUR SUBDOMAIN AND COLLECT DATA
!-----------------------------------------------------------------------
!
          IF(iel.EQ.-2) THEN
            IF(add) THEN
!             A LOST-AGAIN TRACEBACK DETECTED, ALREADY HERE
!             SET THE IMPLANTING PARAMETERS
              iproc=ifapar(ifa  ,elt(iplot))
              iloc =ifapar(ifa+3,elt(iplot))
!             ANOTHER ONE AS IPID, MEANS ALSO NOT LOCALISED
              recvchar(iplot)%NEPID=iproc
              recvchar(iplot)%INE=iloc
            ELSE
              CALL collect_char(ipid,iplot,elt(iplot),ifa,0,0,isp,nsp,
     &                          xplot(iplot),yplot(iplot),0.d0,0.d0,
     &                          dx(iplot),dy(iplot),0.d0,0.d0,
     &                          ifapar,nchdim,nchara)
            ENDIF
!           EXITING LOOP ON ISP
            EXIT
          ENDIF
!
!-----------------------------------------------------------------------
!         SPECIAL TREATMENT FOR SOLID OR LIQUID BOUNDARIES
!-----------------------------------------------------------------------
!
          dxp = dx(iplot)
          dyp = dy(iplot)
          i1  = ikle(elt(iplot),ifa)
          i2  = ikle(elt(iplot),isui(ifa))
          dx1 = x(i2) - x(i1)
          dy1 = y(i2) - y(i1)
!
          IF(iel.EQ.-1) THEN
!
!-----------------------------------------------------------------------
!           HERE SOLID BOUNDARY, VELOCITY IS PROJECTED ON THE BOUNDARY
!           AND WE GO ON
!-----------------------------------------------------------------------
!
!           HERE A1 IS THE PARAMETRIC COORDINATE OF THE PROJECTED
!           DISPLACEMENT ON SEGMENT I1----I2
!
            a1 = (dxp*dx1 + dyp*dy1) / (dx1**2 + dy1**2)
!
!           THE TOTAL DISPLACEMENT IS PROJECTED HERE, NOT THE REMAINING
!           PART, BUT ONLY THE DIRECTION WILL BE USED
            dx(iplot) = a1 * dx1
            dy(iplot) = a1 * dy1
!
!           NOW A1 IS THE PARAMETRIC COORDINATE ON SEGMENT I1----I2
!           OF THE POSITION OF THE ARRIVAL POINT, I.E. INTERSECTION
!           + REMAINING DISPLACEMENT PROJECTED ON THE SEGMENT
!           ITS VALUE MAY BE OUTSIDE THE RANGE (0,1). THE VALUE OF A1
!           SIMPLIFIES INTO THE FOLLOWING FORMULA, BECAUSE IT IS
!           SIMPLY VECTOR I1----P PROJECTED ON SEGMENT I1----I2
!
            a1 = ((xp-x(i1))*dx1+(yp-y(i1))*dy1)/(dx1**2+dy1**2)
            shp(      ifa ,iplot) = 1.d0 - a1
            shp( isui(ifa),iplot) = a1
            shp(isui2(ifa),iplot) = 0.d0
            xplot(iplot) = x(i1) + a1 * dx1
            yplot(iplot) = y(i1) + a1 * dy1
            IF(add) THEN
              recvchar(iplot)%XP=xplot(iplot)
              recvchar(iplot)%YP=yplot(iplot)
              recvchar(iplot)%DX=dx(iplot)
              recvchar(iplot)%DY=dy(iplot)
            ENDIF
!
            GOTO 50
!
          ELSEIF(iel.EQ.0) THEN
!
!------------------------------------------------------------------------
!           HERE WE HAVE A LIQUID BOUNDARY, THE CHARACTERISTIC IS STOPPED
!------------------------------------------------------------------------
!
            denom = dxp*dy1-dyp*dx1
            IF(abs(denom).GT.1.d-12) THEN
              a1  = (dxp*(yp-y(i1))-dyp*(xp-x(i1))) / denom
            ELSE
              a1  = 0.d0
            ENDIF
            a1 = max(min(a1,1.d0),0.d0)
            shp(      ifa ,iplot) = 1.d0 - a1
            shp( isui(ifa),iplot) = a1
            shp(isui2(ifa),iplot) = 0.d0
            xplot(iplot) = x(i1) + a1 * dx1
            yplot(iplot) = y(i1) + a1 * dy1
!           THIS IS A MARKER FOR PARTICLES EXITING A DOMAIN
!           SENS=-1 FOR BACKWARD CHARACTERISTICS
            elt(iplot) = - sens * elt(iplot)
            EXIT
!
          ELSE
!
            WRITE(lu,*) 'UNEXPECTED CASE IN SCHAR11'
            WRITE(lu,*) 'IEL=',iel
            CALL plante(1)
            stop
!
          ENDIF
!
        ENDIF
!
        ENDDO
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE schar11_sto
!                       ******************
                        SUBROUTINE schar12
!                       ******************
!
     &(u,v,dt,nrk,x,y,ikle,ifabor,xplot,yplot,dx,dy,shp,elt,
     & nplot,nelem,nelmax,surdet,sens,
     & ifapar,nchdim,nchara,add)
!
!***********************************************************************
! BIEF VERSION 6.3           24/04/97    J-M JANIN (LNH) 30 87 72 84
!
!***********************************************************************
!
!  FONCTION : THIS IS A MERE COPY OF SCHAR11, EXCEPT THE INTERPOLATION
!             OF VELOCITY WHICH IS HERE CONSIDERED QUASI-BUBBLE
!
!history  J-M HERVOUET (LNHE)
!+        07/04/2013
!+        V6P3
!+   Correct size of velocities given to allow bound checking.
!
!-----------------------------------------------------------------------
!                             ARGUMENTS
! .________________.____.______________________________________________.
! |      NOM       |MODE|                   ROLE                       |
! |________________|____|______________________________________________|
! |    U,V         | -->| COMPOSANTE DE LA VITESSE DU CONVECTEUR       |
! |    DT          | -->| PAS DE TEMPS.                                |
! |    NRK         | -->| NOMBRE DE SOUS-PAS DE RUNGE-KUTTA.           |
! |    X,Y         | -->| COORDONNEES DES POINTS DU MAILLAGE.          |
! |    IKLE        | -->| TRANSITION ENTRE LES NUMEROTATIONS LOCALE    |
! |                |    | ET GLOBALE.                                  |
! |    IFABOR      | -->| NUMEROS DES ELEMENTS AYANT UNE FACE COMMUNE  |
! |                |    | AVEC L'ELEMENT .  SI IFABOR<0 OU NUL         |
! |                |    | ON A UNE FACE LIQUIDE,SOLIDE,OU PERIODIQUE   |
! |  XPLOT,YPLOT   |<-->| POSITIONS SUCCESSIVES DES DERIVANTS.         |
! |    DX,DY       | -- | STOCKAGE DES SOUS-PAS . |
! |    SHP         |<-->| COORDONNEES BARYCENTRIQUES 2D AU PIED DES    |
! |                |    | COURBES CARACTERISTIQUES.                    |
! |    ELT         |<-->| NUMEROS DES ELEMENTS 2D AU PIED DES COURBES  |
! |                |    | CARACTERISTIQUES.                            |
! |    NPLOT       | -->| NOMBRE DE DERIVANTS.                         |
! |    NELEM       | -->| NOMBRE D'ELEMENTS.                           |
! |    NELMAX      | -->| NOMBRE MAXIMAL D'ELEMENTS DANS LE MAILLAGE 2D|
! |    SURDET      | -->| VARIABLE UTILISEE PAR LA TRANSFORMEE ISOPARAM.
! |    SENS        | -->| DESCENTE OU REMONTEE DES CARACTERISTIQUES.   |
! |________________|____|______________________________________________|
!  MODE: -->(DONNEE NON MODIFIEE),<--(RESULTAT),<-->(DONNEE MODIFIEE)
!-----------------------------------------------------------------------
!     - PROGRAMMES APPELES : NEANT
!
!***********************************************************************
!
      USE bief
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)    :: SENS,NCHDIM
      INTEGER         , INTENT(IN)    :: NELEM,NELMAX,NPLOT,NRK
      INTEGER         , INTENT(IN)    :: IKLE(nelmax,4),IFABOR(nelmax,3)
      INTEGER         , INTENT(INOUT) :: ELT(*),NCHARA
      DOUBLE PRECISION, INTENT(IN)    :: U(*),V(*)
      DOUBLE PRECISION, INTENT(IN)    :: SURDET(nelem)
      DOUBLE PRECISION, INTENT(INOUT) :: XPLOT(*),YPLOT(*)
      DOUBLE PRECISION, INTENT(INOUT) :: SHP(3,*)
      DOUBLE PRECISION, INTENT(IN)    :: DT
      DOUBLE PRECISION, INTENT(IN)    :: X(*),Y(*)
      DOUBLE PRECISION, INTENT(INOUT) :: DX(*),DY(*)
      INTEGER, INTENT(IN)             :: IFAPAR(6,*)
      LOGICAL, INTENT(IN)             :: ADD
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER :: IPLOT,ISP,I1,I2,I3,IEL,ISO,IFA,ISUI(3),ISUI2(3)
      INTEGER IPROC,ILOC,ISPDONE,NSP
      DOUBLE PRECISION PAS,A1,DX1,DY1,DXP,DYP,XP,YP,DENOM
      DOUBLE PRECISION SHP11,SHP12,SHP14
      DOUBLE PRECISION SHP22,SHP23,SHP24
      DOUBLE PRECISION SHP33,SHP31,SHP34
!
      INTRINSIC int,max,min,sqrt
!
      parameter( isui = (/ 2 , 3 , 1 /) )
      parameter( isui2 = (/ 3 , 1 , 2 /) )
      DOUBLE PRECISION, PARAMETER :: EPSILO = -1.d-6
!
!-----------------------------------------------------------------------
!     FOR EVERY POINT
!-----------------------------------------------------------------------
!
      DO iplot=1,nplot
!
        IF(add) THEN
!
          xplot(iplot)   = recvchar(iplot)%XP
          yplot(iplot)   = recvchar(iplot)%YP
          dx(iplot)      = recvchar(iplot)%DX
          dy(iplot)      = recvchar(iplot)%DY
          elt(iplot)     = recvchar(iplot)%INE
          nsp            = recvchar(iplot)%NSP ! R-K STEPS TO BE FULLFILLED
          ispdone        = recvchar(iplot)%ISP ! R-K STEPS ALREADY DONE
          iel = elt(iplot)
          xp  = xplot(iplot)
          yp  = yplot(iplot)
          i1 = ikle(iel,1)
          i2 = ikle(iel,2)
          i3 = ikle(iel,3)
          shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &           -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
          shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &           -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
          shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &           -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
!         ASSUME TO BE LOCALISED IT WILL BE SET OTHERWISE IF LOST-AGAIN
          recvchar(iplot)%NEPID=-1
!
        ELSE
          iel = elt(iplot)
!         POINTS WITH IEL=0 ARE TREATED SO THAT THE FINAL INTERPOLATION
!         GIVES 0., AND WE SKIP TO NEXT POINT IPLOT (CYCLE)
!         THIS WILL NOT INTERFERE WITH ELT(IPLOT)=0 GIVEN ON LIQUID BOUNDARIES
!         BY ARRAY IFABOR, THAT MAY HAPPEN LATER
          IF(iel.EQ.0) THEN
            elt(iplot)=1
            shp(1,iplot)=0.d0
            shp(2,iplot)=0.d0
            shp(3,iplot)=0.d0
            cycle
          ENDIF
          shp11=shp(1,iplot)-shp(3,iplot)
          shp12=shp(2,iplot)-shp(3,iplot)
          shp14=3.d0*shp(3,iplot)
          shp22=shp(2,iplot)-shp(1,iplot)
          shp23=shp(3,iplot)-shp(1,iplot)
          shp24=3.d0*shp(1,iplot)
          shp33=shp(3,iplot)-shp(2,iplot)
          shp31=shp(1,iplot)-shp(2,iplot)
          shp34=3.d0*shp(2,iplot)
          IF(     shp11.GT.     2.d0*epsilo .AND.
     &            shp11.LT.1.d0-4.d0*epsilo .AND.
     &            shp12.GT.     2.d0*epsilo .AND.
     &            shp12.LT.1.d0-4.d0*epsilo .AND.
     &            shp14.LT.1.d0-4.d0*epsilo ) THEN
            dxp = u(ikle(iel,1)) * shp11
     &          + u(ikle(iel,2)) * shp12
     &          + u(ikle(iel,4)) * shp14
            dyp = v(ikle(iel,1)) * shp11
     &          + v(ikle(iel,2)) * shp12
     &          + v(ikle(iel,4)) * shp14
          ELSEIF( shp22.GT.     2.d0*epsilo .AND.
     &            shp22.LT.1.d0-4.d0*epsilo .AND.
     &            shp23.GT.     2.d0*epsilo .AND.
     &            shp23.LT.1.d0-4.d0*epsilo .AND.
     &            shp24.LT.1.d0-4.d0*epsilo ) THEN
            dxp = u(ikle(iel,2)) * shp22
     &          + u(ikle(iel,3)) * shp23
     &          + u(ikle(iel,4)) * shp24
            dyp = v(ikle(iel,2)) * shp22
     &          + v(ikle(iel,3)) * shp23
     &          + v(ikle(iel,4)) * shp24
          ELSEIF( shp33.GT.     2.d0*epsilo .AND.
     &            shp33.LT.1.d0-4.d0*epsilo .AND.
     &            shp31.GT.     2.d0*epsilo .AND.
     &            shp31.LT.1.d0-4.d0*epsilo .AND.
     &            shp34.LT.1.d0-4.d0*epsilo ) THEN
            dxp = u(ikle(iel,3)) * shp33
     &          + u(ikle(iel,1)) * shp31
     &          + u(ikle(iel,4)) * shp34
            dyp = v(ikle(iel,3)) * shp33
     &          + v(ikle(iel,1)) * shp31
     &          + v(ikle(iel,4)) * shp34
          ELSE
            WRITE(lu,*) 'SCHAR12: POINT ',iplot
            WRITE(lu,*) '         NOT IN ELEMENT ',iel
            WRITE(lu,*) 'SHP(1,IPLOT)=',shp(1,iplot)
            WRITE(lu,*) 'SHP(2,IPLOT)=',shp(2,iplot)
            WRITE(lu,*) 'SHP(3,IPLOT)=',shp(3,iplot)
            WRITE(lu,*) 'EPSILO=',epsilo,' IPID=',ipid,' CASE 1'
            WRITE(lu,*) 'SHP11=',shp11
            WRITE(lu,*) 'SHP12=',shp12
            WRITE(lu,*) 'SHP14=',shp14
            WRITE(lu,*) 'SHP22=',shp22
            WRITE(lu,*) 'SHP23=',shp23
            WRITE(lu,*) 'SHP24=',shp24
            WRITE(lu,*) 'SHP33=',shp33
            WRITE(lu,*) 'SHP31=',shp31
            WRITE(lu,*) 'SHP34=',shp34
            CALL plante(1)
            stop
          ENDIF
          nsp=max(1,int(nrk*dt*sqrt((dxp**2+dyp**2)*surdet(iel))))
          ispdone=1
        ENDIF
!
        pas = sens*dt/nsp
!
!       LOOP ON RUNGE-KUTTA SUB-STEPS
!
!       COMPILER MUST DO NOTHING IF ISPDONE>NSP
!       IN MODE "ADD", ISP = ISPDONE HAS NOT BEEN FULLY DONE
!       IT IS RESTARTED HERE
!
        DO isp=ispdone,nsp
!
!-----------------------------------------------------------------------
!       LOCALISING THE ARRIVAL POINT
!-----------------------------------------------------------------------
!
!       IN MODE "ADD" ITERATIONS ALREADY DONE ARE SKIPPED AND
!                     CHARACTERISTICS GONE IN ANOTHER SUB-DOMAIN SKIPPED
!
        IF(add) THEN
          IF(isp.EQ.ispdone) GO TO 50
          IF(recvchar(iplot)%NEPID.NE.-1) cycle
        ENDIF
!
        iel = elt(iplot)
        i1 = ikle(iel,1)
        i2 = ikle(iel,2)
        i3 = ikle(iel,3)
!
        shp11=shp(1,iplot)-shp(3,iplot)
        shp12=shp(2,iplot)-shp(3,iplot)
        shp14=3.d0*shp(3,iplot)
        shp22=shp(2,iplot)-shp(1,iplot)
        shp23=shp(3,iplot)-shp(1,iplot)
        shp24=3.d0*shp(1,iplot)
        shp33=shp(3,iplot)-shp(2,iplot)
        shp31=shp(1,iplot)-shp(2,iplot)
        shp34=3.d0*shp(2,iplot)
        IF(     shp11.GT.     2.d0*epsilo .AND.
     &          shp11.LT.1.d0-4.d0*epsilo .AND.
     &          shp12.GT.     2.d0*epsilo .AND.
     &          shp12.LT.1.d0-4.d0*epsilo .AND.
     &          shp14.LT.1.d0-4.d0*epsilo ) THEN
          dx(iplot) = ( u(ikle(iel,1)) * shp11
     &                + u(ikle(iel,2)) * shp12
     &                + u(ikle(iel,4)) * shp14 ) * pas
          dy(iplot) = ( v(ikle(iel,1)) * shp11
     &                + v(ikle(iel,2)) * shp12
     &                + v(ikle(iel,4)) * shp14 ) * pas
        ELSEIF( shp22.GT.     2.d0*epsilo .AND.
     &          shp22.LT.1.d0-4.d0*epsilo .AND.
     &          shp23.GT.     2.d0*epsilo .AND.
     &          shp23.LT.1.d0-4.d0*epsilo .AND.
     &          shp24.LT.1.d0-4.d0*epsilo ) THEN
          dx(iplot) = ( u(ikle(iel,2)) * shp22
     &                + u(ikle(iel,3)) * shp23
     &                + u(ikle(iel,4)) * shp24 ) * pas
          dy(iplot) = ( v(ikle(iel,2)) * shp22
     &                + v(ikle(iel,3)) * shp23
     &                + v(ikle(iel,4)) * shp24 ) * pas
        ELSEIF( shp33.GT.     2.d0*epsilo .AND.
     &          shp33.LT.1.d0-4.d0*epsilo .AND.
     &          shp31.GT.     2.d0*epsilo .AND.
     &          shp31.LT.1.d0-4.d0*epsilo .AND.
     &          shp34.LT.1.d0-4.d0*epsilo ) THEN
          dx(iplot) = ( u(ikle(iel,3)) * shp33
     &                + u(ikle(iel,1)) * shp31
     &                + u(ikle(iel,4)) * shp34 ) * pas
          dy(iplot) = ( v(ikle(iel,3)) * shp33
     &                + v(ikle(iel,1)) * shp31
     &                + v(ikle(iel,4)) * shp34 ) * pas
        ELSE
          WRITE(lu,*) 'SCHAR12: POINT ',iplot
          WRITE(lu,*) '         NOT IN ELEMENT ',elt(iplot)
          WRITE(lu,*) 'SHP(1,IPLOT)=',shp(1,iplot)
          WRITE(lu,*) 'SHP(2,IPLOT)=',shp(2,iplot)
          WRITE(lu,*) 'SHP(3,IPLOT)=',shp(3,iplot)
          WRITE(lu,*) 'EPSILO=',epsilo,' IPID=',ipid,' CASE 2'
          CALL plante(1)
          stop
        ENDIF
        xp = xplot(iplot) + dx(iplot)
        yp = yplot(iplot) + dy(iplot)
!
        shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                 -(y(i3)-y(i2))*(xp-x(i2))) * surdet(iel)
        shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                 -(y(i1)-y(i3))*(xp-x(i3))) * surdet(iel)
        shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                 -(y(i2)-y(i1))*(xp-x(i1))) * surdet(iel)
!
        xplot(iplot) = xp
        yplot(iplot) = yp
!
        IF(add) THEN
!         CONTINUOUS SETTING OF THE REACHED POSITION FOR IPLOT
!         AND THE NUMBER OF STEPS DONE ALREADY
          recvchar(iplot)%XP=xplot(iplot)
          recvchar(iplot)%YP=yplot(iplot)
          recvchar(iplot)%DX=dx(iplot)
          recvchar(iplot)%DY=dy(iplot)
          recvchar(iplot)%INE=elt(iplot)
          recvchar(iplot)%ISP=isp
        ENDIF
!
!-----------------------------------------------------------------------
!       TEST: IS THE PATHLINE WENT OUT THE ORIGINAL ELEMENT
!-----------------------------------------------------------------------
!
50      CONTINUE
!
        iso = 0
        IF(shp(1,iplot).LT.epsilo) iso = 1
        IF(shp(2,iplot).LT.epsilo) iso = iso + 2
        IF(shp(3,iplot).LT.epsilo) iso = iso + 4
!
        IF(iso.NE.0) THEN
!
!-----------------------------------------------------------------------
!         HERE WE ARE OUT OF THE ELEMENT
!-----------------------------------------------------------------------
!
          iel = elt(iplot)
          xp = xplot(iplot)
          yp = yplot(iplot)
!
          IF(iso.EQ.1) THEN
            ifa = 2
          ELSEIF (iso.EQ.2) THEN
            ifa = 3
          ELSEIF (iso.EQ.4) THEN
            ifa = 1
          ELSEIF (iso.EQ.3) THEN
            ifa = 2
            IF(dx(iplot)*(y(ikle(iel,3))-yp).LT.
     &                dy(iplot)*(x(ikle(iel,3))-xp)) ifa = 3
            ELSEIF (iso.EQ.6) THEN
              ifa = 3
              IF (dx(iplot)*(y(ikle(iel,1))-yp).LT.
     &            dy(iplot)*(x(ikle(iel,1))-xp)) ifa = 1
            ELSE
              ifa = 1
              IF(dx(iplot)*(y(ikle(iel,2))-yp).LT.
     &           dy(iplot)*(x(ikle(iel,2))-xp)) ifa = 2
            ENDIF
!
            iel = ifabor(iel,ifa)
!
            IF(iel.GT.0) THEN
!
!-----------------------------------------------------------------------
!             HERE WE ARRIVE IN ANOTHER ELEMENT
!-----------------------------------------------------------------------
!
              i1 = ikle(iel,1)
              i2 = ikle(iel,2)
              i3 = ikle(iel,3)
!
              elt(iplot) = iel
              shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &                       -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
              shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &                       -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
              shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &                       -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
!
              GOTO 50
!
            ENDIF
!
!-----------------------------------------------------------------------
!           HERE WE PASS TO NEIGHBOUR SUBDOMAIN AND COLLECT DATA
!-----------------------------------------------------------------------
!
            IF(iel.EQ.-2) THEN
              IF(add) THEN
!               A LOST-AGAIN TRACEBACK DETECTED, ALREADY HERE
!               SET THE IMPLANTING PARAMETERS
                iproc=ifapar(ifa  ,elt(iplot))
                iloc =ifapar(ifa+3,elt(iplot))
!               ANOTHER ONE AS IPID, MEANS ALSO NOT LOCALISED
                recvchar(iplot)%NEPID=iproc
                recvchar(iplot)%INE=iloc
              ELSE
                CALL collect_char(ipid,iplot,elt(iplot),ifa,0,0,isp,nsp,
     &                            xplot(iplot),yplot(iplot),0.d0,0.d0,
     &                            dx(iplot),dy(iplot),0.d0,0.d0,
     &                            ifapar,nchdim,nchara)
              ENDIF
!             EXITING LOOP ON ISP
              EXIT
            ENDIF
!
!-----------------------------------------------------------------------
!           SPECIAL TREATMENT FOR SOLID OR LIQUID BOUNDARIES
!-----------------------------------------------------------------------
!
            dxp = dx(iplot)
            dyp = dy(iplot)
            i1  = ikle(elt(iplot),ifa)
            i2  = ikle(elt(iplot),isui(ifa))
            dx1 = x(i2) - x(i1)
            dy1 = y(i2) - y(i1)
!
            IF(iel.EQ.-1) THEN
!
!-----------------------------------------------------------------------
!             HERE SOLID BOUNDARY, VELOCITY IS PROJECTED ON THE BOUNDARY
!             AND WE GO ON
!-----------------------------------------------------------------------
!
              a1 = (dxp*dx1 + dyp*dy1) / (dx1**2 + dy1**2)
              dx(iplot) = a1 * dx1
              dy(iplot) = a1 * dy1
!
              a1 = ((xp-x(i1))*dx1+(yp-y(i1))*dy1)/(dx1**2+dy1**2)
              shp(      ifa ,iplot) = 1.d0 - a1
              shp( isui(ifa),iplot) = a1
              shp(isui2(ifa),iplot) = 0.d0
              xplot(iplot) = x(i1) + a1 * dx1
              yplot(iplot) = y(i1) + a1 * dy1
              IF(add) THEN
                recvchar(iplot)%XP=xplot(iplot)
                recvchar(iplot)%YP=yplot(iplot)
                recvchar(iplot)%DX=dx(iplot)
                recvchar(iplot)%DY=dy(iplot)
              ENDIF
!
              GOTO 50
!
            ELSEIF(iel.EQ.0) THEN
!
!-----------------------------------------------------------------------
!             HERE WE HAVE A LIQUID BOUNDARY, THE CHARACTERISTIC IS STOPPED
!-----------------------------------------------------------------------
!
              denom = dxp*dy1-dyp*dx1
              IF(abs(denom).GT.1.d-12) THEN
                a1  = (dxp*(yp-y(i1))-dyp*(xp-x(i1))) / denom
              ELSE
                a1  = 0.d0
              ENDIF
              a1 = max(min(a1,1.d0),0.d0)
              shp(      ifa ,iplot) = 1.d0 - a1
              shp( isui(ifa),iplot) = a1
              shp(isui2(ifa),iplot) = 0.d0
              xplot(iplot) = x(i1) + a1 * dx1
              yplot(iplot) = y(i1) + a1 * dy1
!             THIS IS A MARKER FOR PARTICLES EXITING A DOMAIN
!             SENS=-1 FOR BACKWARD CHARACTERISTICS
              elt(iplot) = - sens * elt(iplot)
              EXIT
!
            ELSE
!
              WRITE(lu,*) 'UNEXPECTED CASE IN SCHAR12'
              WRITE(lu,*) 'IEL=',iel
              CALL plante(1)
              stop
!
            ENDIF
!
          ENDIF
!
        ENDDO
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE schar12
!                       ******************
                        SUBROUTINE schar13
!                       ******************
!
     &(u,v,dt,nrk,x,y,ikle,ifabor,xplot,yplot,dx,dy,shp,elt,
     & nplot,nelem,nelmax,surdet,sens,
     & ifapar,nchdim,nchara,add)
!
!***********************************************************************
! BIEF VERSION 6.2           24/04/97    J-M JANIN (LNH) 30 87 72 84
!
!***********************************************************************
!
!  FONCTION : THIS IS A MERE COPY OF SCHAR11, EXCEPT THE INTERPOLATION
!             OF VELOCITY WHICH IS HERE CONSIDERED QUADRATIC
!
!history  J-M HERVOUET (LNHE)
!+        07/04/2013
!+        V6P3
!+   Size of velocities set to * to blind bound checking. The correct
!+   quadratic size would be needed here.
!
!-----------------------------------------------------------------------
!                             ARGUMENTS
! .________________.____.______________________________________________.
! |      NOM       |MODE|                   ROLE                       |
! |________________|____|______________________________________________|
! |    U,V         | -->| COMPOSANTE DE LA VITESSE DU CONVECTEUR       |
! |    DT          | -->| PAS DE TEMPS.                                |
! |    NRK         | -->| NOMBRE DE SOUS-PAS DE RUNGE-KUTTA.           |
! |    X,Y         | -->| COORDONNEES DES POINTS DU MAILLAGE.          |
! |    IKLE        | -->| TRANSITION ENTRE LES NUMEROTATIONS LOCALE    |
! |                |    | ET GLOBALE.                                  |
! |    IFABOR      | -->| NUMEROS DES ELEMENTS AYANT UNE FACE COMMUNE  |
! |                |    | AVEC L'ELEMENT .  SI IFABOR<0 OU NUL         |
! |                |    | ON A UNE FACE LIQUIDE,SOLIDE,OU PERIODIQUE   |
! |  XPLOT,YPLOT   |<-->| POSITIONS SUCCESSIVES DES DERIVANTS.         |
! |    DX,DY       | -- | STOCKAGE DES SOUS-PAS . |
! |    SHP         |<-->| COORDONNEES BARYCENTRIQUES 2D AU PIED DES    |
! |                |    | COURBES CARACTERISTIQUES.                    |
! |    ELT         |<-->| NUMEROS DES ELEMENTS 2D AU PIED DES COURBES  |
! |                |    | CARACTERISTIQUES.                            |
! |    NPLOT       | -->| NOMBRE DE DERIVANTS.                         |
! |    NPOIN       | -->| NOMBRE DE POINTS DU MAILLAGE.                |
! |    NELEM       | -->| NOMBRE D'ELEMENTS.                           |
! |    NELMAX      | -->| NOMBRE MAXIMAL D'ELEMENTS DANS LE MAILLAGE 2D|
! |    SURDET      | -->| VARIABLE UTILISEE PAR LA TRANSFORMEE ISOPARAM.
! |    SENS        | -->| DESCENTE OU REMONTEE DES CARACTERISTIQUES.   |
! |________________|____|______________________________________________|
!  MODE: -->(DONNEE NON MODIFIEE),<--(RESULTAT),<-->(DONNEE MODIFIEE)
!-----------------------------------------------------------------------
!     - PROGRAMMES APPELES : NEANT
!
!***********************************************************************
!
      USE bief
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)    :: SENS,NCHDIM
      INTEGER         , INTENT(IN)    :: NELEM,NELMAX,NPLOT,NRK
      INTEGER         , INTENT(IN)    :: IKLE(nelmax,6),IFABOR(nelmax,3)
      INTEGER         , INTENT(INOUT) :: ELT(nplot),NCHARA
!                                        QUADRATIC VELOCITIES
      DOUBLE PRECISION, INTENT(IN)    :: U(*),V(*),SURDET(nelem)
      DOUBLE PRECISION, INTENT(INOUT) :: XPLOT(nplot),YPLOT(nplot)
      DOUBLE PRECISION, INTENT(INOUT) :: SHP(3,nplot)
      DOUBLE PRECISION, INTENT(IN)    :: DT
      DOUBLE PRECISION, INTENT(IN)    :: X(*),Y(*)
      DOUBLE PRECISION, INTENT(INOUT) :: DX(nplot),DY(nplot)
      INTEGER, INTENT(IN)             :: IFAPAR(6,*)
      LOGICAL, INTENT(IN)             :: ADD
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IPLOT,ISP,I1,I2,I3,I4,I5,I6,IEL,ISO,IFA
      INTEGER :: ISUI(3),ISUI2(3)
      INTEGER IPROC,ILOC,ISPDONE,NSP
      DOUBLE PRECISION PAS,A1,DX1,DY1,DXP,DYP,XP,YP,DENOM
!
      INTRINSIC int,max,min,sqrt
!
      parameter( isui =  (/ 2 , 3 , 1 /) )
      parameter( isui2 = (/ 3 , 1 , 2 /) )
      DOUBLE PRECISION, PARAMETER :: EPSILO = -1.d-6
!
!-----------------------------------------------------------------------
!     FOR EVERY POINT
!-----------------------------------------------------------------------
!
      DO iplot=1,nplot
!
        IF(add) THEN
!
          xplot(iplot)   = recvchar(iplot)%XP
          yplot(iplot)   = recvchar(iplot)%YP
          dx(iplot)      = recvchar(iplot)%DX
          dy(iplot)      = recvchar(iplot)%DY
          elt(iplot)     = recvchar(iplot)%INE
          nsp            = recvchar(iplot)%NSP ! R-K STEPS TO BE FULLFILLED
          ispdone        = recvchar(iplot)%ISP ! R-K STEPS ALREADY DONE
          iel = elt(iplot)
          xp  = xplot(iplot)
          yp  = yplot(iplot)
          i1 = ikle(iel,1)
          i2 = ikle(iel,2)
          i3 = ikle(iel,3)
          shp(1,iplot) = ((x(i3)-x(i2))*(yp-y(i2))
     &           -(y(i3)-y(i2))*(xp-x(i2)))*surdet(iel)
          shp(2,iplot) = ((x(i1)-x(i3))*(yp-y(i3))
     &           -(y(i1)-y(i3))*(xp-x(i3)))*surdet(iel)
          shp(3,iplot) = ((x(i2)-x(i1))*(yp-y(i1))
     &           -(y(i2)-y(i1))*(xp-x(i1)))*surdet(iel)
!         ASSUME TO BE LOCALISED IT WILL BE SET OTHERWISE IF LOST-AGAIN
          recvchar(iplot)%NEPID=-1
!
        ELSE
          iel = elt(iplot)
!         POINTS WITH IEL=0 ARE TREATED SO THAT THE FINAL INTERPOLATION
!         GIVES 0., AND WE SKIP TO NEXT POINT IPLOT (CYCLE)