derlag.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE derlag
!                   *****************
!
     &(u,v,dt,x,y,lt,ielm,ielmu,ndp,npoin,nelem,nelmax,xlag,ylag,dx,dy,
     & nsp,shplag,deblag,finlag,eltlag,nlag,resux,resuy,ispdone,mesh)
!
!***********************************************************************
! BIEF   V6P2                                   21/08/2010
!***********************************************************************
!
!brief    - SETS THE BARYCENTRIC COORDINATES IN THE MESH,
!+                  AT THE START OF COMPUTATION FOR EACH DRIFTING FLOAT.
!+                  HERE WE COMPUTE THE LAGRANGIAN DRIFT.
!+
!+            - COMPUTES THE SUCCESSIVE POSITIONS OF THIS FLOAT
!+                 (SUBSEQUENT TIMESTEPS).
!
!warning  Will not work in parallel (this would require calling scaract
!+        instead of char11, and adaptation of scaract)
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  J-M HERVOUET (LNHE)
!+        19/06/2012
!+        V6P2
!+   Adapted to call SCARACT instead of CHAR11. However further
!+   modifications are required for parallelism.
!
!history  J-M HERVOUET (LNHE)
!+        28/10/2014
!+        V7P0
!+   Stopping drifts that exit the domain has been added
!+   + error message : parallelism does not work here.
!+   (possible but a long implementation copied on particles)
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        14/08/2015
!+        V7P1
!+   Hardcoded argument NRK added for SCARACT.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| DEBLAG         |-->| TIME STEP FOR STARTING THE COMPUTATION
!| DT             |-->| TIME STEP
!| DX             |<->| WORK ARRAY
!| DY             |<->| WORK ARRAY
!| ELTLAG         |<->| ELEMENT NUMBERS OF FLOATS
!| FINLAG         |-->| TIME STEP FOR ENDING THE COMPUTATION
!| IELM           |-->| TYPE OF ELEMENT IN THE MESH
!| IELMU          |-->| TYPE OF ELEMENT FOR THE VELOCITIES
!| LT             |-->| TIME STEP NUMBER.
!| MASKEL         |-->| MASKING OF ELEMENTS.
!|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
!| MASKPT         |-->| MASKING PER POINT.
!| MESH           |-->| MESH STRUCTURE
!| NDP            |-->| NUMBER OF POINTS PER ELEMENT
!| NELEM          |-->| NUMBER OF ELEMENTS.
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS.
!| NLAG           |-->| NOMBER OF FLOATS.
!| NPOIN          |-->| NUMBER OF POINTS
!| NSP            |-->| NUMBER OF SUB-STEPS IN THE RUNGE-KUTTA METHOD
!| RESUX          |<--| ARRAY WITH SUCCESSIVE ABSCISSAE OF FLOATS
!| RESUY          |<--| ARRAY WITH SUCCESSIVE ORDINATES OF FLOATS
!| SHPLAG         |<->| BARYCENTRIC COORDINATES OF FLOATS
!|                |   | IN THEIR ELEMENTS.
!| T8             |-->| BLOCK OF WORK BIEF_OBJ STRUCTURES.
!| U              |-->| X-COMPONENT OF VELOCITY
!| V              |-->| Y-COMPONENT OF VELOCITY
!| X              |-->| ABSCISSAE OF POINTS IN THE MESH
!| XLAG           |<->| INSTANTANEOUS X POSITIONS OF FLOATS
!| Y              |-->| ORDINATES OF POINTS IN THE MESH
!| YLAG           |<->| INSTANTANEOUS Y POSITIONS OF FLOATS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_derlag => derlag
      USE streamline
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)    :: NPOIN,LT,IELM,NDP,NELEM,NLAG
      INTEGER         , INTENT(IN)    :: NELMAX,IELMU
      DOUBLE PRECISION, INTENT(IN)    :: U(npoin),V(npoin),DT
      DOUBLE PRECISION, INTENT(IN)    :: X(npoin),Y(npoin)
      DOUBLE PRECISION, INTENT(INOUT) :: XLAG(npoin,nlag)
      DOUBLE PRECISION, INTENT(INOUT) :: YLAG(npoin,nlag)
      INTEGER         , INTENT(INOUT) :: DEBLAG(nlag),FINLAG(nlag)
      INTEGER         , INTENT(INOUT) :: ELTLAG(npoin,nlag)
      INTEGER         , INTENT(INOUT) :: ISPDONE(npoin)
      DOUBLE PRECISION, INTENT(INOUT) :: DX(npoin),DY(npoin)
      INTEGER         , INTENT(INOUT) :: NSP(npoin)
      DOUBLE PRECISION, INTENT(INOUT) :: RESUX(npoin),RESUY(npoin)
      DOUBLE PRECISION, INTENT(INOUT) :: SHPLAG(ndp,npoin,nlag)
      TYPE(bief_mesh) , INTENT(INOUT) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER ILAG,JLAG,LTT,IPOIN,ETA(1),SENS,NPLAN,FRE(1),FREBUF(1),NRK
      TYPE(bief_obj) :: SVOID
!
      DOUBLE PRECISION ZSTAR(1),ZCONV(1),SHZ(1),Z(1),SHF(1)
      DOUBLE PRECISION SHPBUF(3,1),SHZBUF(1)
!
!-----------------------------------------------------------------------
!
      IF(ncsize.GT.1) THEN
        WRITE(lu,*) 'LAGRANGIAN DRIFTS NOT PROGRAMMED IN PARALLEL'
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
!     HARDCODED NUMBER OF SUB-STEPS FOR COMPUTING THE PATH-LINES
!
      nrk=3
!
!     2D HERE
!
      nplan=1
!
!     FORWARD CHARACTERISTICS
!
      sens=1
!
      DO ilag=1,nlag
!
        IF(lt.EQ.deblag(ilag)) THEN
!
!-----------------------------------------------------------------------
!
!   - SETS THE BARYCENTRIC COORDINATES IN THE MESH , AT THE START
!     OF COMPUTATION FOR EACH FLOAT
!
!-----------------------------------------------------------------------
!
          IF(ielm.EQ.11) THEN
!
!  P1 TRIANGLES
!  ============
!
!      FILLS THE SHP AND ELT (OPTIMISED)
!
            CALL gtsh11(shplag(1,1,ilag),eltlag(1,ilag),mesh%IKLE%I,
     &                  mesh%ELTCAR%I,
     &                  npoin,nelem,nelmax,1,.false.,.false.)
!                                          1=NSEG, WRONG VALUE, NOT USED
!
          ELSE
!
            WRITE(lu,*) ielm,': ELEMENT NOT IMPLEMENTED IN DERLAG'
            CALL plante(1)
            stop
!
          ENDIF
!
          CALL ov('X=Y     ', x=xlag(1,ilag), y=x, dim1=npoin)
          CALL ov('X=Y     ', x=ylag(1,ilag), y=y, dim1=npoin)
!
        ELSEIF(lt.GT.deblag(ilag).AND.lt.LE.finlag(ilag)) THEN
!
!-----------------------------------------------------------------------
!
!     COMPUTES THE SUCCESSIVE POSITIONS OF THIS FLOAT
!     (SUBSEQUENT TIMESTEPS)
!
!-----------------------------------------------------------------------
!
          CALL scaract(svoid,svoid,u,v,v,v,x,y,
     &                 zstar,zstar,xlag(1,ilag),ylag(1,ilag),
     &                 zconv,zconv,
     &                 dx,dy,dy,dy,z,shplag(1,1,ilag),shz,shf,
     &                 mesh%SURDET%R,
     &                 dt,mesh%IKLE%I,mesh%IFABOR%I,eltlag(1,ilag),
     &                 eta,fre,nsp,ispdone,ielm,ielmu,nelem,nelmax,
     &                 0,npoin,ndp,nrk,nplan,1,
     &                 mesh,npoin,bief_nbpts(ielmu,mesh),sens,
!                      PROVISIONAL, THIS WILL NOT WORK IN PARALLEL
     &                 shpbuf,shzbuf,shzbuf,frebuf,1,
     &                 .true.)
!                       POST : DATA KEPT FOR A POSTERIORI INTERPOLATION
!                              HERE FOR NEXT STEP
!
!-----------------------------------------------------------------------
!
!       STOPPING THE FLOATS THAT LEFT THE DOMAIN
!
!-----------------------------------------------------------------------
!
          DO ipoin=1,npoin
            IF(eltlag(ipoin,ilag).LT.0) THEN
!             THIS POINT WILL NO LONGER BE TREATED
              eltlag(ipoin,ilag) = 0
!             NEXT TWO LINES IN PREVIOUS VERSIONS
!             BUT WHY NOT KEEPING THE LAST POSITION AT THE EXIT?
!             XLAG(IPOIN,ILAG) = X(IPOIN)
!             YLAG(IPOIN,ILAG) = Y(IPOIN)
            ENDIF
          ENDDO
!
!-----------------------------------------------------------------------
!
        ENDIF
!
      ENDDO ! ILAG
!
!-----------------------------------------------------------------------
!
!   - STORAGE FOR RESULTS OUTPUT OF THE LAST COMPUTED FLOAT
!
!-----------------------------------------------------------------------
!
      CALL ov('X=C     ', x=resux, c=0.d0, dim1=npoin)
      CALL ov('X=C     ', x=resuy, c=0.d0, dim1=npoin)
      ltt=0
      jlag=1
      DO ilag=1,nlag
        IF(finlag(ilag).GT.ltt.AND.finlag(ilag).LE.lt) THEN
          ltt=finlag(ilag)
          jlag=ilag
        ENDIF
      ENDDO
      IF(ltt.NE.0) THEN
        CALL ov('X=Y-Z   ', x=resux, y=xlag(1,jlag), z=x, dim1=npoin)
        CALL ov('X=Y-Z   ', x=resuy, y=ylag(1,jlag), z=y, dim1=npoin)
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END