topogr.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE topogr
!                   *****************
!
     &(zf,zref,zfe,ikle,ifabor,nbor,nelbor,nulone,
     & itra05,itra02,itra03,nelem,nptfr,npoin,mxptvs)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    FINELY ANALYSES THE TOPOGRAPHY AND BUILDS ZFE.
!+
!+            THE ARRAY OF BOTTOM ELEVATIONS BY ELEMENTS: ZFE
!+                WILL ENSURE IN THE FUTURE THAT THERE WILL NOT BE
!+                LIQUID DOMAINS CONNECTED BY A SINGLE NODE.
!
!history  J-M JANIN (LNH)
!+        17/08/94
!+        V5P1
!+
!
!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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IFABOR         |-->| ELEMENTS BEHIND THE EDGES OF A TRIANGLE
!|                |   | IF NEGATIVE OR ZERO, THE EDGE IS A LIQUID
!|                |   | BOUNDARY
!| IKLE           |-->| CONNECTIVITY TABLE.
!| ITRA02         |<--| INTEGER WORK ARRAY
!| ITRA03         |<--| INTEGER WORK ARRAY
!| ITRA05         |<--| INTEGER WORK ARRAY
!| MXPTVS         |-->| MAXIMUM NUMBER OF NEIGHBOURS OF A POINT
!| NBOR           |-->| GLOBAL NUMBER OF BOUNDARY POINTS
!| NELBOR         |-->| FOR THE KTH BOUNDARY EDGE, GIVES THE CORRESPONDING
!|                |   | ELEMENT.
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NPOIN          |-->| NUMBER OF POINTS
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| NULONE         |-->| GOES WITH ARRAY NELBOR. NELBOR GIVES THE
!|                |   | ADJACENT ELEMENT, NULONE GIVES THE LOCAL
!|                |   | NUMBER OF THE FIRST NODE OF THE BOUNDARY EDGE
!|                |   | I.E. 1, 2 OR 3 FOR TRIANGLES.
!| ZF             |-->| ELEVATION OF BOTTOM, PER POINT
!| ZFE            |<--| ELEVATION OF BOTTOM, PER ELEMENT
!| ZREF           |<--| CORRECTED BOTTOM ELEVATION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)    :: NELEM,NPTFR,NPOIN,MXPTVS
      INTEGER, INTENT(IN)    :: IKLE(nelem,3),IFABOR(nelem,3)
      INTEGER, INTENT(IN)    :: NBOR(nptfr),NELBOR(nptfr),NULONE(nptfr)
      INTEGER, INTENT(INOUT) :: ITRA05(npoin),ITRA02(npoin)
      INTEGER, INTENT(INOUT) :: ITRA03(npoin)
      DOUBLE PRECISION, INTENT(INOUT) :: ZFE(nelem),ZREF(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: ZF(npoin)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELEM,IPTFR,IPOIN,I,I1,I2,I3,N1,N2,ERR,IMAX
      LOGICAL FLAG
!
!     DYNAMICALLY ALLOCATES INTEGERS
!
      INTEGER, DIMENSION(:,:), ALLOCATABLE :: ITRA01,ITRA04,IFAN
!
      INTEGER :: IPREV(3)
      parameter( iprev = (/ 3 , 1 , 2 /) )
      DOUBLE PRECISION, PARAMETER :: EPSILO = 1.d-6
!
!-----------------------------------------------------------------------
!
      ALLOCATE(ifan(nelem,3)          ,stat=err)
      CALL check_allocate(err, 'IFAN')
      ALLOCATE(itra01(npoin,mxptvs+1) ,stat=err)
      CALL check_allocate(err, 'ITRA01')
      ALLOCATE(itra04(npoin,6)        ,stat=err)
      CALL check_allocate(err, 'ITRA04')
!
!-----------------------------------------------------------------------
!
!    FILLS IN IFAN
!
!   IFAN(IELEM,IFACE) GIVES THE LOCAL NUMBER IN THE ELEMENT DEFINED
!   BY IFABOR(IELEM,IFACE) OF THE NODE WITH LOCAL NUMBER IFACE IN
!   ELEMENT IELEM.
!
!    FIRST GO AT FILLING ZFE
!
!    STARTS TO FILL IN ITRA01 AND ITRA02
!
!   ITRA01 AND ITRA02 PERFORM THE REVERSE OPERATION FROM IKLE.
!   IN THIS LOOP ITRA01(IPOIN,1) IS FILLED IN. IT GIVES THE BIGGEST
!   ELEMENT NUMBER CONTAINING IPOIN. ITRA02(IPOIN) GIVES THE LOCAL
!   NUMBER OF IPOIN IN THIS ELEMENT.
!
!-----------------------------------------------------------------------
!
!  INITIALISES ITRA01 TO DETECT HOLES IN THE MESH
!  (POINTS NOT LINKED TO AN ELEMENT, CASE OF CURVILINEAR MESH)
!
      DO i1    = 1 , mxptvs+1
        DO ipoin = 1 , npoin
          itra01(ipoin,i1) = 0
        ENDDO ! IPOIN
      ENDDO ! I1
!
!
      DO ielem = 1,nelem
!
        i1 = ikle(ielem,1)
        i2 = ikle(ielem,2)
        i3 = ikle(ielem,3)
!
        ifan(ielem,1) = 0
        n1 = ifabor(ielem,1)
        IF (n1.GT.0) THEN
          IF(ikle(n1,1).EQ.i1) ifan(ielem,1) = 3
          IF(ikle(n1,2).EQ.i1) ifan(ielem,1) = 1
          IF(ikle(n1,3).EQ.i1) ifan(ielem,1) = 2
        ENDIF
!
        ifan(ielem,2) = 0
        n1 = ifabor(ielem,2)
        IF (n1.GT.0) THEN
          IF(ikle(n1,1).EQ.i2) ifan(ielem,2) = 3
          IF(ikle(n1,2).EQ.i2) ifan(ielem,2) = 1
          IF(ikle(n1,3).EQ.i2) ifan(ielem,2) = 2
        ENDIF
!
        ifan(ielem,3) = 0
        n1 = ifabor(ielem,3)
        IF (n1.GT.0) THEN
          IF(ikle(n1,1).EQ.i3) ifan(ielem,3) = 3
          IF(ikle(n1,2).EQ.i3) ifan(ielem,3) = 1
          IF(ikle(n1,3).EQ.i3) ifan(ielem,3) = 2
        ENDIF
!
        zfe(ielem) = max(zf(i1),zf(i2),zf(i3))
        itra01(i1,1) = ielem
        itra02(i1)   = 1
        itra01(i2,1) = ielem
        itra02(i2)   = 2
        itra01(i3,1) = ielem
        itra02(i3)   = 3
!
      ENDDO ! IELEM
!
!-----------------------------------------------------------------------
!
!    STARTS TO FILL IN ITRA01 AND ITRA02 FOR BOUNDARY NODES
!
!   FOR THESE POINTS ITRA01(IPOIN,1) IS NOT ANY ELEMENT: IT'S THE
!   ELEMENT CONTAINING A BOUNDARY SIDE BETWEEN NODE IPOIN AND THE
!   FOLLOWING NODE IN THE LOCAL NUMBERING OF THE ELEMENT.
!
!-----------------------------------------------------------------------
!
      DO iptfr = 1,nptfr
        itra01(nbor(iptfr),1) = nelbor(iptfr)
        itra02(nbor(iptfr))   = nulone(iptfr)
      ENDDO ! IPTFR
!
!-----------------------------------------------------------------------
!
!    RESUMES TO FILL IN ITRA01
!
!   ITRA01(IPOIN,I+1) IS THE NUMBER OF THE ELEMENT ADJACENT TO ELEMENT
!   ITRA01(IPOIN,1) WHEN TURNING ANTI-CLOCKWISE AROUND POINT IPOIN.
!   ITRA02 IS A UNIDIMENSIONAL ARRAY BECAUSE THE LOCAL NUMBER OF IPOIN
!   IN ANY ELEMENT ITRA01(IPOIN,*) WILL NOT NEEDED IN THE FUTURE.
!
!
!    FILLS IN ITRA03
!
!   ITRA03(IPOIN) CORRESPONDS TO THE NUMBER OF ELEMENTS CONTAINING IPOIN
!   (WITH A NEGATIVE SIGN IF IPOIN IS A BOUNDARY NODE).
!
!    BEWARE |||
!
!   ONE POINT CAN ONLY BE PART OF A MAXIMUM OF 10 ELEMENTS
!
!-----------------------------------------------------------------------
!
      DO ipoin = 1,npoin
        itra03(ipoin) = 0
      ENDDO ! IPOIN
!
      imax = 0
!
40    CONTINUE
      flag = .false.
      imax = imax + 1
      IF (imax.GT.mxptvs+1) THEN
        WRITE(lu,24) mxptvs
24      FORMAT(1x,'TOPOGR : THE MAXIMUM NUMBER OF NEIGHBOURS TO'/,1x,
     &            '         A POINT IS GREATER THAN THE VALUE  ',/,1x,
     &            '         GIVEN BY MXPTVS :',1i6)
        CALL plante(0)
        stop
      ENDIF
!
      DO ipoin = 1,npoin
!
        IF (itra03(ipoin).EQ.0) THEN
          n1 = itra01(ipoin,imax)
          IF(n1.NE.0) THEN
            flag = .true.
            n2 = ifabor(n1,iprev(itra02(ipoin)))
!                        HERE IMAX IS NEVER AT ITS MAXIMUM
            itra01(ipoin,imax+1) = n2
            itra02(ipoin) = ifan(n1,iprev(itra02(ipoin)))
            IF (n2.LE.0)               itra03(ipoin) = -imax
            IF (n2.EQ.itra01(ipoin,1)) itra03(ipoin) =  imax
          ENDIF
        ENDIF
!
      ENDDO ! IPOIN
!
      IF (flag) GOTO 40
!
60    CONTINUE
!
!-----------------------------------------------------------------------
!
!    DETERMINES LOCAL EXTREMA FOR ZFE BY TURNING AROUND A NODE
!
!   ITRA04(IPOIN,I) CORRESPONDS TO THE IEME EXTREMUM FOUND, WITH THE
!   ASSOCIATED ELEMENT GIVEN BY ITRA01(IPOIN,ITRA04(IPOIN,I)).
!   ITRA02(IPOIN) GIVES THE TOTAL NUMBER OF INCREASE AND DECREASE STAGES
!   (WITH A NEGATIVE SIGN IF THE LAST STAGE FOUND IS A DECREASE).
!
!
!-----------------------------------------------------------------------
!
      DO ipoin = 1,npoin
        itra02(ipoin) = 0
        itra05(ipoin) = 0
      ENDDO ! IPOIN
!
      DO i = 1,imax-1
!
        DO ipoin = 1,npoin
!
          IF (itra03(ipoin).GE.i.OR.itra03(ipoin).LT.-i) THEN
!
            n1 = itra01(ipoin,i)
            n2 = itra01(ipoin,i+1)
!
            IF (zfe(n2).GT.zfe(n1)+epsilo) THEN
              IF (itra02(ipoin).LT.0) itra04(ipoin,-itra02(ipoin))=i
              IF (itra02(ipoin).LE.0) itra02(ipoin)=-itra02(ipoin)+1
            ELSEIF (zfe(n2).LT.zfe(n1)-epsilo) THEN
              IF (itra02(ipoin).GT.0) itra04(ipoin, itra02(ipoin))=i
              IF (itra02(ipoin).GE.0) itra02(ipoin)=-itra02(ipoin)-1
            ENDIF
!
          ENDIF
!
        ENDDO ! IPOIN
!
      ENDDO ! I
!
      DO ipoin = 1,npoin
        IF((itra03(ipoin).LT.0.AND.(itra02(ipoin).LE.-4.OR.
     &      itra02(ipoin).GE.5)).OR.abs(itra02(ipoin)).GE.6) THEN
          WRITE(lu,*) 'THE MESH AROUND THE NODE ',ipoin,' HAS TO'
          WRITE(lu,*) 'BE REFINED BECAUSE OF THE BATHYMETRY'
          CALL plante(1)
          stop
        ENDIF
      ENDDO ! IPOIN
!
!-----------------------------------------------------------------------
!
!    CORRECTS ZFE DEPENDING ON ITRA02
!
!-----------------------------------------------------------------------
!
      flag = .false.
!
      DO ipoin = 1,npoin
!
        i1 = itra03(ipoin)
!
        IF (i1.LT.0) THEN
!
!-----------------------------------------------------------------------
!
!    CORRECTS ZFE FOR BOUNDARY NODES
!
!   IF ITRA02(IPOIN) EQUALS 1, -1, 2 : NO CORRECTION
!   IF ITRA02(IPOIN) EQUALS -2, 3, -3, 4 : CORRECTION
!
!-----------------------------------------------------------------------
!
          IF (itra02(ipoin).EQ.-2) THEN
!
            flag = .true.
            IF (zfe(itra01(ipoin,-i1)).GT.zfe(itra01(ipoin,1))) THEN
              itra02(ipoin) = itra04(ipoin,1) + 1
              itra05(ipoin) = -i1
            ELSE
              itra02(ipoin) = 1
              itra05(ipoin) = itra04(ipoin,1) - 1
            ENDIF
            zref(ipoin) = zfe(itra01(ipoin,itra04(ipoin,1)))
!
          ELSEIF (itra02(ipoin).EQ.3) THEN
!
            flag = .true.
            IF (zfe(itra01(ipoin,itra04(ipoin,2))).GT.
     &          zfe(itra01(ipoin,1))) THEN
              itra02(ipoin) = itra04(ipoin,1) + 1
              itra05(ipoin) = -i1
            ELSE
              itra02(ipoin) = 1
              itra05(ipoin) = itra04(ipoin,1) - 1
            ENDIF
            zref(ipoin) = zfe(itra01(ipoin,itra04(ipoin,1)))
!
          ELSEIF (itra02(ipoin).EQ.-3) THEN
!
            flag = .true.
            IF (zfe(itra01(ipoin,-i1)).GT.
     &          zfe(itra01(ipoin,itra04(ipoin,1)))) THEN
              itra02(ipoin) = itra04(ipoin,2) + 1
              itra05(ipoin) = -i1
            ELSE
              itra02(ipoin) = 1
              itra05(ipoin) = itra04(ipoin,2) - 1
            ENDIF
            zref(ipoin) = zfe(itra01(ipoin,itra04(ipoin,2)))
!
          ELSEIF (itra02(ipoin).EQ.4) THEN
!
            flag = .true.
            IF (zfe(itra01(ipoin,itra04(ipoin,3))).GT.
     &          zfe(itra01(ipoin,itra04(ipoin,1)))) THEN
              itra02(ipoin) = itra04(ipoin,2) + 1
              itra05(ipoin) = -i1
            ELSE
              itra02(ipoin) = 1
              itra05(ipoin) = itra04(ipoin,2) - 1
            ENDIF
            zref(ipoin) = zfe(itra01(ipoin,itra04(ipoin,2)))
!
          ENDIF
!
        ELSE
!
!-----------------------------------------------------------------------
!
!    CORRECTS ZFE FOR INTERIOR NODES
!
!   IF ITRA02(IPOIN) EQUALS 1, -1, 2, -2, 3, -3 : NO CORRECTION
!   IF ITRA02(IPOIN) EQUALS 4, -4, 5, -5 : CORRECTION
!
!-----------------------------------------------------------------------
!
          IF (itra02(ipoin).EQ.4) THEN
!
            flag = .true.
            IF (zfe(itra01(ipoin,itra04(ipoin,3))).GT.
     &          zfe(itra01(ipoin,itra04(ipoin,1)))) THEN
              itra02(ipoin) = itra04(ipoin,2) + 1
              itra05(ipoin) = i1
            ELSE
              itra02(ipoin) = 2
              itra05(ipoin) = itra04(ipoin,2) - 1
            ENDIF
            zref(ipoin) = min(zfe(itra01(ipoin,itra04(ipoin,2))),
     &                        zfe(itra01(ipoin,1)))
!
          ELSEIF (itra02(ipoin).EQ.-4) THEN
!
            flag = .true.
            IF (zfe(itra01(ipoin,itra04(ipoin,2))).GT.
     &          zfe(itra01(ipoin,1))) THEN
              itra02(ipoin) = itra04(ipoin,1) + 1
              itra05(ipoin) = itra04(ipoin,3) - 1
            ELSE
              itra02(ipoin) = mod(itra04(ipoin,3),i1) + 1
              itra05(ipoin) = itra04(ipoin,1) - 1
            ENDIF
            zref(ipoin) = min(zfe(itra01(ipoin,itra04(ipoin,1))),
     &                        zfe(itra01(ipoin,itra04(ipoin,3))))
!
          ELSEIF (itra02(ipoin).EQ.5) THEN
!
            flag = .true.
            IF (zfe(itra01(ipoin,itra04(ipoin,4))).GT.
     &          zfe(itra01(ipoin,itra04(ipoin,2)))) THEN
              itra02(ipoin) = itra04(ipoin,3) + 1
              itra05(ipoin) = itra04(ipoin,1) - 1
            ELSE
              itra02(ipoin) = itra04(ipoin,1) + 1
              itra05(ipoin) = itra04(ipoin,3) - 1
            ENDIF
            zref(ipoin) = min(zfe(itra01(ipoin,itra04(ipoin,1))),
     &                        zfe(itra01(ipoin,itra04(ipoin,3))))
!
          ELSEIF (itra02(ipoin).EQ.-5) THEN
!
            flag = .true.
            IF (zfe(itra01(ipoin,itra04(ipoin,3))).GT.
     &          zfe(itra01(ipoin,itra04(ipoin,1)))) THEN
              itra02(ipoin) = itra04(ipoin,2) + 1
              itra05(ipoin) = itra04(ipoin,4) - 1
            ELSE
              itra02(ipoin) = mod(itra04(ipoin,4),i1) + 1
              itra05(ipoin) = itra04(ipoin,2) - 1
            ENDIF
            zref(ipoin) = min(zfe(itra01(ipoin,itra04(ipoin,2))),
     &                        zfe(itra01(ipoin,itra04(ipoin,4))))
!
          ENDIF
!
        ENDIF
!
      ENDDO ! IPOIN
!
      IF (flag) THEN
!
        DO ipoin = 1,npoin
!
          IF (itra05(ipoin).NE.0) THEN
!
            IF (itra05(ipoin).LT.itra02(ipoin)) THEN
              DO i = itra02(ipoin),itra03(ipoin)
                zfe(itra01(ipoin,i)) = max(zfe(itra01(ipoin,i)),
     &                                     zref(ipoin))
              ENDDO ! I
              itra02(ipoin) = 1
            ENDIF
!
            DO i = itra02(ipoin),itra05(ipoin)
              zfe(itra01(ipoin,i)) = max(zfe(itra01(ipoin,i)),
     &                                   zref(ipoin))
            ENDDO ! I
!
          ENDIF
!
        ENDDO ! IPOIN
!
        GOTO 60
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      DEALLOCATE(ifan)
      DEALLOCATE(itra01)
      DEALLOCATE(itra04)
!
!-----------------------------------------------------------------------
!
      RETURN
      END