read_drogues.f

Go to the documentation of this file.

!                   ***********************
                    SUBROUTINE read_drogues
!                   ***********************
!
     &( np,np_max,xp,yp,zp,tagp,clsp,eltp,shpp,
     &  npoin2,npoin3,nelem,nelmax,ikle,x,y,t2dplo )
!
!***********************************************************************
! BIEF   V8P0
!***********************************************************************
!
!brief    Reads the results for a given time step from a drogues file.
!+        The drogues file will have to have been created previously
!+        by the TELEMAC system
!
!history  S.E. BOURBAN (HRW)
!+        21/08/2018
!+        V8P0
!+        New subroutine created to allow for continuous simulation.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CLSP           |<--| CLASS OF DROGUES
!| ELTP           |<--| ELEMENT FOR EACH DROGUE
!| IKLE           |-->| CONNECTIVITY TABLE
!| NELEM          |-->| NUMBER OF ELEMENTS IN THE MESH
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS IN THE MESH
!| NP             |<--| NUMBER OF DROGUES
!| NPOIN2         |-->| NUMBER OF POINTS OF 2D MESH
!| NPOIN3         |-->| NUMBER OF POINTS OF 3D MESH
!| NP_MAX         |-->| MAXIMUM NUMBER OF DROGUES
!| SHPP           |<--| BARYCENTRIC COORDINATES OF DROGUES IN THEIR
!|                |   | ELEMENTS.
!| TAGP           |<--| TAG OF DROGUES
!| T2DPLO         |-->| PREVIOUS DROGUES COMPUTATION FILE
!| X              |-->| ABSCISSAE OF POINTS IN THE MESH
!| XP             |<--| ABSCISSAE OF DROGUES
!| Y              |-->| ORDINATES OF POINTS IN THE MESH
!| YP             |<--| ORDINATES OF DROGUES
!| ZP             |<--| ELEVATIONS OF DROGUES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE initial_drogues, ONLY : ndrg_tags
      USE interface_parallel, ONLY : p_max
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_file) , INTENT(IN) :: T2DPLO
      INTEGER         , INTENT(IN)    :: NP_MAX,NELEM,NELMAX
      INTEGER         , INTENT(IN)    :: NPOIN2,NPOIN3
      INTEGER         , INTENT(INOUT) :: NP
      INTEGER, INTENT(IN)             :: IKLE(nelmax,3)
      INTEGER, INTENT(INOUT)          :: TAGP(np_max),ELTP(np_max)
      INTEGER, INTENT(INOUT)          :: CLSP(np_max)
      DOUBLE PRECISION, INTENT(IN)    :: X(npoin2),Y(npoin2)
      DOUBLE PRECISION, INTENT(INOUT) :: XP(np_max),YP(np_max)
      DOUBLE PRECISION, INTENT(INOUT) :: SHPP(3,np_max),ZP(np_max)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, PARAMETER      :: SIZELIGN = 3000
      CHARACTER(LEN=SIZELIGN) :: LIGNE
!
      INTEGER          LUPLO, I,J,ILINE,ICHAR, ITIME(9), IFLOT,MFLOT
      INTEGER           JPID, N1,N2,N3
      LOGICAL          FOUND
      DOUBLE PRECISION  X0,Y0,X1,Y1,X2,Y2,X3,Y3,DET1,DET2,DET3,SURDET
      REAL              :: RVALUE
      REAL, ALLOCATABLE :: RARRAY(:)
      DOUBLE PRECISION, PARAMETER :: CHOUIA = 1.d-10
!
!-----------------------------------------------------------------------
!
!     HEADER OF BLUE KENUE PARCEL FILE IN ASCII FORM
!
      IF( t2dplo%FMT .EQ.'BKBINPCL' .AND. t2dplo%LU.NE.0 ) THEN
        luplo = t2dplo%LU
        CLOSE(luplo)
!
!-----------------------------------------------------------------------
!
!     > READING THE HEADER OF THE DROGUES FILE
!
        OPEN( luplo, file=trim(t2dplo%TELNAME),form='FORMATTED',
     &        access='STREAM', action='READ' )
        iline = 0
        ichar = 0
        found = .false.
        rewind(luplo)
!       COUNTS COMMENTS FOR NOW (TODO: RECORD ATTRIBUTES)
        DO WHILE( .true. )
          iline = iline + 1
          READ( luplo,fmt='(A)',end=1001,err=1901 ) ligne
          ichar = ichar + len(trim(ligne)) + 2
          IF( ligne(1:18).EQ.':WrittenBy TELEMAC' ) found = .true.
          IF( ligne(1:10).EQ.':EndHeader' ) EXIT
        ENDDO
        CLOSE(luplo)
!
!       CORE OF BLUE KENUE PARCEL FILE IN BINARY LITTLE ENDIAN FORM
        OPEN( luplo, file=trim(t2dplo%TELNAME),form='UNFORMATTED',
     &        action='READ',convert='LITTLE_ENDIAN',
     &        access='STREAM', status='OLD' )
!
!       LOOKING FOR THE MAXIMUM NUMBER OF DROGUES FOR MEMORY ALLOCATION
        mflot = 0
        iline = 0
        rewind(luplo)
        READ(luplo) ligne(1:ichar)
        DO WHILE( .true. )
          READ(luplo,end=23,err=1902) itime
          iline = iline + 1
          READ(luplo,err=1902) iflot
          mflot = max( mflot,iflot )
          READ(luplo) ( rvalue, i = 1,iflot*4 )
          IF( npoin2.NE.npoin3 ) THEN
            READ(luplo) ( rvalue, i = 1,iflot )
          ENDIF
          READ(luplo) iflot
        ENDDO
 23     CONTINUE
!
!       CHECK AGAINST NUMBER OF PARCELS POSSIBLE
        IF( mflot.GT.np_max ) THEN
          WRITE(lu,33) mflot,np_max
 33       FORMAT(1x,'READ_DROGUES:',/,
     &    1x,'     REQUIRED NUMBER OF DROGUES (',i8,')',/,
     &    1x,'     LARGER THAN THE MAXIMUM NUMBER OF DROGUES',/,
     &    1x,'     POSSIBLE (',i8,')',/,
     &    1x,'     INCREASE THE MAXIMUM OR REDUCE YOUR DENSITY.')
          CALL plante(1)
          stop
        ENDIF
!
!-----------------------------------------------------------------------
!
!     > READING ALL RECORDS BUT THE LAST ONE
!
        ALLOCATE( rarray(mflot) )
        rewind(luplo)
!
        READ(luplo) ligne(1:ichar)
        DO j = 1,iline-1
          READ(luplo) itime
          READ(luplo) iflot
          READ(luplo) ( rvalue, i = 1,iflot*4 )
          IF( npoin2.NE.npoin3 ) THEN
            READ(luplo) ( rvalue, i = 1,iflot )
          ENDIF
          READ(luplo) iflot
        ENDDO
!
!-----------------------------------------------------------------------
!
!     > READING THE LAST RECORD
!
        READ(luplo) itime
        READ(luplo) iflot
!       X
        READ(luplo) ( rarray(i), i = 1,iflot )
        DO i = 1,iflot
          xp(i) = rarray(i)
        ENDDO
!       Y
        READ(luplo) ( rarray(i), i = 1,iflot )
        DO i = 1,iflot
          yp(i) = rarray(i)
        ENDDO
!       Z
        IF( npoin2.NE.npoin3 ) THEN
          READ(luplo) ( rarray(i), i = 1,iflot )
          DO i = 1,iflot
            zp(i) = rarray(i)
          ENDDO
        ENDIF
!       TAG
        READ(luplo) ( rarray(i), i = 1,iflot )
        DO i = 1,iflot
          tagp(i) = int( rarray(i) )
        ENDDO
!       CLASS
        READ(luplo) ( rarray(i), i = 1,iflot )
        DO i = 1,iflot
          clsp(i) = int( rarray(i) )
        ENDDO
!
        DEALLOCATE(rarray)
        CLOSE(luplo)
!
!-----------------------------------------------------------------------
!
!     > POSITIONING THE DORGUES WITHIN THE MESH
!
        np = 0
        IF( iflot.GT.0 ) THEN
          DO i = 1,iflot
!
            x0 = xp(i)
            y0 = yp(i)
!
!       IDENTIFYING THE TRIANGLE WHERE THE PARCELS FELL IN
!       (FOR EACH PROCESSOR, AND IT MAY BE ON THE EDGE OF TWO OR MORE)
!
            jpid = -1
            DO j = 1,nelem
!
              n1 = ikle( j,1 )
              n2 = ikle( j,2 )
              n3 = ikle( j,3 )
              x1 = x(n1)
              y1 = y(n1)
              x2 = x(n2)
              y2 = y(n2)
              x3 = x(n3)
              y3 = y(n3)
!
              det1 = ( x3-x2 )*( y0-y2 ) - ( y3-y2 )*( x0-x2 )
              det2 = ( x1-x3 )*( y0-y3 ) - ( y1-y3 )*( x0-x3 )
              det3 = ( x2-x1 )*( y0-y1 ) - ( y2-y1 )*( x0-x1 )
              IF( det1.GE.-chouia .AND.
     &            det2.GE.-chouia .AND.
     &            det3.GE.-chouia ) THEN
                jpid = ipid
                EXIT
              ENDIF
!
            ENDDO
!
!        /!\ IT IS IMPORTANT TO NOTE THAT NP .LE. IFLOT
!
            IF( ncsize.GT.1 ) THEN
              jpid = p_max(jpid)
            ENDIF
            IF( jpid.EQ.ipid ) THEN
              np = np + 1
              xp(np) = x0
              yp(np) = y0
              eltp(np) = j
              surdet = 1.d0 / ( (x2-x1)*(y3-y1) - (x3-x1)*(y2-y1) )
              shpp(1,np) = det1*surdet
              shpp(2,np) = det2*surdet
              shpp(3,np) = det3*surdet
              tagp(np) = tagp(i)
              clsp(np) = clsp(i)
            ENDIF
!
          ENDDO
        ENDIF
!
        ndrg_tags = np
        IF( ncsize.GT.1 ) ndrg_tags = p_max(ndrg_tags)
        WRITE(lu,34) np
 34     FORMAT(1x,'READ_DROGUES:',/,i8,
     &    1x,'     DROGUES HAVE BEEN READ FROM THE',/,
     &    1x,'     PREVIOUS DROGUES FILE.')
!
!-----------------------------------------------------------------------
!
      ENDIF
      RETURN
!
!-----------------------------------------------------------------------
!
 1001 CONTINUE
      WRITE(lu,*) ligne
      WRITE(lu,*) 'REACHED THE END OF THE DROGUES FILE'
      WRITE(lu,*) 'BEFORE FINDING THE END OF THE HEADER.'
      WRITE(lu,*) 'MAYBE THIS IS NOT A DROGUES FILE.'
      CALL plante(1)
      stop
 1901 CONTINUE
      WRITE(lu,*) 'READING ERROR ON THE DROGUES FILE'
      WRITE(lu,*) 'AT LINE OF DATA : ',iline
      CALL plante(1)
      stop
 1902 CONTINUE
      WRITE(lu,*) 'READING ERROR ON THE DROGUES FILE'
      WRITE(lu,*) 'AFTER RECORD : ',itime
      CALL plante(1)
      stop
!
!-----------------------------------------------------------------------
!
      RETURN
      END