charac.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE charac
!                   *****************
!
     &( fn    , ftild  , nomb  , uconv  , vconv , wconv  , frconv ,
     &  zstar , freq   ,
     &  dt    , ifamas , ielm  , npoin2 , nplan , jf     , nf     ,
     &  msk   , shp   , shz    , shf   , tb     , elt    ,
     &  eta   , fre    , it3   , isub   , frebuf, mesh   ,
     &  nelem2, nelmax2, ikle2 , surdet2, am1   , rhs    , slv    ,
     &  agglo , listin , ngauss, unsv   , optcha, post   , perio  ,
     &  ya4d  , sigma  , stocha, visc )
!
!***********************************************************************
! BIEF   V6P3                                   21/08/2010
!***********************************************************************
!
!brief    CALLS THE METHOD OF CHARACTERISTICS
!+               (SUBROUTINE CARACT).
!
!warning  So far the size of some arrays cannot be checked:
!+        FRE, ELT, ETA, IT3, ISUB, FREBUF, it would be better to send
!+        BIEF_OBJ structures.
!
!history  J-M HERVOUET (LNHE)
!+        12/02/2010
!+        V6P0
!+   First version.
!
!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)
!+        03/01/2012
!+        V6P2
!+   NPOIN instead of NPOIN2 in the call to SCARACT at the position
!+   of argument NPLOT (goes with corrections in Streamline.f)
!
!history  J-M HERVOUET (LNHE)
!+        18/04/2013
!+        V6P3
!+   NPOIN instead of NPOIN2 in the call to SCARACT at the position
!+   of argument NPLOT (goes with corrections in Streamline.f)
!
!history  J-M HERVOUET (LNHE)
!+        02/05/2013
!+        V6P3
!+   Arguments STOCHA and VISC added. Attribute target added for some
!+   others (it helps prepare weak form of characteristics).
!
!history  J-M HERVOUET (LNHE)
!+        16/05/2013
!+        V6P3
!+   Attribute target for ISUB.
!
!history  J-M HERVOUET (LNHE)
!+        18/03/2013
!+        V6P3
!+   Now weak form of characteristics also available, see OPTCHA.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AGGLO          |-->| MASS-LUMPING (FOR WEAK FORM)
!| AM1            |<->| A MATRIX (FOR WEAK FORM)
!| RHS            |<->| A RIGHT-HAND SIDE (FOR WEAK FORM)
!| DT             |-->| TIME STEP
!| ELT            |<->| ARRIVAL ELEMENT
!| ETA            |<->| ARRIVAL LAYER (IN 3D WITH PRISMS)
!| FRCONV         |-->| FREQUENCY COMPONENT OF ADVECTION FIELD
!| FRE            |<->| ARRIVAL FREQUENCY (IN 4D)
!| FREBUF         |<->| INTEGER WORK ARRAY (IN 4D)
!| FREQ           |-->| DISCRETISED FREQUENCIES (IN 4D).
!|                |   | IF NOT TOMAWAC, MUST BE ZSTAR !!!!!!!!!!!
!| FN             |-->| VARIABLES AT TIME N .
!| FTILD          |<--| VARIABLES AFTER ADVECTION .
!| IELM           |-->| TYPE OF ELEMENT : 11 : TRIANGLE P1
!|                |   |                   41 : PRISM IN TELEMAC3D
!| IFAMAS         |-->| A MODIFIED IFABOR WHEN ELEMENTS ARE MASKED
!| IKLE2          |-->| CONNECTIVITY TABLE FOR TRIANGLES
!| IT3            |<->| INTEGER WORK ARRAY
!| ISUB           |<->| ARRIVAL SUB-DOMAIN (IN PARALLEL)
!| JF             |-->| FREQUENCY (IN A RANGE OF 1 TO NF)
!| LISTIN         |-->| IF YES, PRINTS INFORMATIONS ON LISTING (WEAK FORM)
!| MESH           |-->| MESH STRUCTURE
!| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
!| NELEM2         |-->| NUMBER OF ELEMENTS IN 2D
!| NELMAX2        |-->| MAXIMUM NUMBER OF ELEMENTS IN 2D
!| NF             |-->| NUMBER OF FREQUENCIES (IN 4D)
!| NGAUSS         |-->| NUMBER OF GAUSS POINTS (WEAK FORM)
!| NOMB           |-->| NUMBER OF VARIABLES TO BE ADVECTED
!| NPLAN          |-->| NUMBER OF PLANES IN THE 3D MESH OF PRISMS
!| NPOIN2         |-->| NUMBER OF POINTS IN THE 2D MESH
!| OPTCHA         |-->| OPTION FOR THE FORM OF CHARACTERISTICS
!|                |   | 1 : STRONG (CLASSICAL) FORM
!|                |   | 2 : BACKWARD WEAK FORM
!| PERIO          |-->| IF YES, PERIODIC VERSION ON THE VERTICAL
!| POST           |-->| IF YES, DATA MUST BE KEPT FOR A POSTERIORI
!|                |   | INTERPOLATION
!| SHP            |<->| BARYCENTRIC COORDINATES OF POINTS IN TRIANGLES
!| SHZ            |<->| BARYCENTRIC COORDINATES ON VERTICAL
!| SHF            |<->| BARYCENTRIC COORDINATES ON THE FREQUENCY AXIS
!| SIGMA          |-->| IF YES, TRANSFORMES MESH FOR TELEMAC-3D
!| SLV            |-->| A SOLVER CONFIGURATION (FOR WEAK FORM)
!| STOCHA         |-->| STOCHASTIC DIFFUSION MODEL
!|                |   | 0: NO DIFFUSION 1: BROWNIAN MODEL
!| SURDET2        |-->| GEOMETRIC COEFFICIENT USED IN PARAMETRIC TRANSFORMATION
!| TB             |<->| BLOCK CONTAINING THE BIEF_OBJ WORK ARRAYS
!| UCONV          |-->| X-COMPONENT OF ADVECTION FIELD
!| UNSV           |-->| 1/(INTEGRAL OF TEST FUNCTIONS)
!| VCONV          |-->| Y-COMPONENT OF ADVECTION FIELD
!| VISC           |-->| VISCOSITY (MAY BE TENSORIAL)
!| WCONV          |-->| Z-COMPONENT OF ADVECTION FIELD IN THE TRANSFORMED MESH
!| YA4D           |-->| IF YES, 4D VERSION FOR TOMAWAC
!| ZSTAR          |-->| TRANSFORMED VERTICAL COORDINATES IN 3D
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_charac => charac
      USE streamline, ONLY : scaract
      USE declarations_telemac, ONLY: deja_charac,t1weak,t2weak,
     &                                t3weak,t4weak,t5weak,t6weak,
     &                                t7weak,shpwea,ftild_weak,shpbuf,
     &                                shzbuf,shzwea
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)           :: NOMB,OPTCHA,NGAUSS
      INTEGER         , INTENT(IN)           :: NPLAN,JF,NF,NELEM2
      INTEGER         , INTENT(IN)           :: NPOIN2,NELMAX2
      INTEGER         , INTENT(INOUT)        :: IELM,FRE(*)
!     NEXT 4 DIMENSIONS ARE A MINIMUM IT IS MORE WITH WEAK FORM
      INTEGER         , INTENT(INOUT)        :: ELT(*)
      INTEGER         , INTENT(INOUT)        :: ETA(*)
      INTEGER         , INTENT(INOUT)        :: IT3(*)
      INTEGER         , INTENT(INOUT)        :: ISUB(*)
      INTEGER         , INTENT(INOUT)        :: FREBUF(*)
      TYPE(bief_obj)  , INTENT(IN)           :: FN,UCONV,VCONV,WCONV
      TYPE(bief_obj)  , INTENT(IN)           :: FRCONV
      TYPE(bief_obj)  , INTENT(IN)           :: ZSTAR,IKLE2
      TYPE(bief_obj)  , INTENT(IN)           :: SURDET2,FREQ,UNSV
      TYPE(bief_obj)  , INTENT(INOUT)        :: TB,SHF,AM1,RHS
      TYPE(bief_obj)  , INTENT(INOUT),TARGET :: FTILD,SHP,SHZ
      LOGICAL         , INTENT(IN)           :: MSK,LISTIN
      DOUBLE PRECISION, INTENT(IN)           :: DT,AGGLO
      TYPE(bief_mesh) , INTENT(INOUT)        :: MESH
      TYPE(bief_obj)  , INTENT(IN), TARGET   :: IFAMAS
      TYPE(slvcfg)    , INTENT(INOUT)        :: SLV
!
      LOGICAL, INTENT(IN), OPTIONAL          :: POST,PERIO,YA4D,SIGMA
      INTEGER, INTENT(IN), OPTIONAL          :: STOCHA
      TYPE(bief_obj), INTENT(IN), OPTIONAL, TARGET :: VISC
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER NPOIN,IELMU,SIZEBUF,ASTOCHA
!
!-----------------------------------------------------------------------
!
      TYPE(bief_obj), POINTER :: T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,AVISC
      TYPE(bief_obj), POINTER :: PT_FTILD,PT_SHPBUF
      INTEGER, DIMENSION(:), POINTER :: IFA
      DOUBLE PRECISION, DIMENSION(:), POINTER :: PT_SHP,PT_SHZ
      INTEGER I,NPT,DIM1F,IPLAN,NG,NRK
      LOGICAL QUAD,QUAB,APOST,APERIO,AYA4D,ASIGMA
      INTRINSIC min
!
!-----------------------------------------------------------------------
!
!     HARDCODED NUMBER OF SUB-STEPS FOR COMPUTING THE PATH-LINES
!
      nrk=3
!
!-----------------------------------------------------------------------
!
      IF(optcha.GT.1) THEN
        IF(ielm.EQ.11) THEN
          ng=ngauss*nelem2
        ELSEIF(ielm.EQ.41) THEN
          ng=ngauss*nelem2*(nplan-1)
        ENDIF
        IF(.NOT.deja_charac) THEN
          CALL allblo(ftild_weak,'FTIWEA')
          CALL bief_allvec_in_block(ftild_weak,ftild%N,1,
     &                              'FTW   ',ng,1,0,mesh)
          CALL bief_allvec(1,t1weak,'T1WEAK',ng,1,0,mesh)
          CALL bief_allvec(1,t2weak,'T2WEAK',ng,1,0,mesh)
          CALL bief_allvec(1,t4weak,'T4WEAK',ng,1,0,mesh)
          CALL bief_allvec(1,t5weak,'T5WEAK',ng,1,0,mesh)
          CALL bief_allvec(1,shpwea,'SHPWEA',ng,3,0,mesh)
          IF(ielm.EQ.41) THEN
            CALL bief_allvec(1,t3weak,'T3WEAK',ng,1,0,mesh)
            CALL bief_allvec(1,t6weak,'T6WEAK',ng,1,0,mesh)
            CALL bief_allvec(1,t7weak,'T7WEAK',ng,1,0,mesh)
            CALL bief_allvec(1,shzwea,'SHZWEA',ng,1,0,mesh)
          ELSE
            CALL bief_allvec(1,t3weak,'T3WEAK', 1,1,0,mesh)
            CALL bief_allvec(1,t6weak,'T6WEAK', 1,1,0,mesh)
            CALL bief_allvec(1,t7weak,'T7WEAK', 1,1,0,mesh)
            CALL bief_allvec(1,shzwea,'SHZWEA',1 ,1,0,mesh)
          ENDIF
          IF(ncsize.GT.1) THEN
            CALL bief_allvec(1,shpbuf,'SHPBUF',ng,3,0,mesh)
          ELSE
            CALL bief_allvec(1,shpbuf,'SHPBUF',1 ,3,0,mesh)
          ENDIF
          IF(ncsize.GT.1.AND.ielm.EQ.41) THEN
            CALL bief_allvec(1,shzbuf,'SHZBUF',ng,1,0,mesh)
          ELSE
            CALL bief_allvec(1,shzbuf,'SHZBUF',1 ,1,0,mesh)
          ENDIF
          deja_charac=.true.
        ENDIF
      ENDIF
!
!-----------------------------------------------------------------------
!  OPTIONAL OPTIONS
!-----------------------------------------------------------------------
!
!     ENABLING A POSTERIORI INTERPOLATION
!
      IF(PRESENT(post)) THEN
        apost=post
      ELSE
        apost=.false.
      ENDIF
!
!     PERIODICITY FOR TOMAWAC
!
      IF(PRESENT(perio)) THEN
        aperio=perio
      ELSE
        aperio=.false.
      ENDIF
!
!     4D FOR TOMAWAC
!
      IF(PRESENT(ya4d)) THEN
        aya4d=ya4d
      ELSE
        aya4d=.false.
      ENDIF
!
!     TRANSFORMED MESH FOR TELEMAC-3D
!
      IF(PRESENT(sigma)) THEN
        asigma=sigma
      ELSE
        asigma=.false.
      ENDIF
!
!     STOCHASTIC DIFFUSION
!
      IF(PRESENT(stocha)) THEN
        astocha=stocha
      ELSE
        astocha=0
      ENDIF
      IF(PRESENT(visc)) THEN
        avisc => visc
      ELSE
        IF(astocha.NE.0) THEN
          WRITE(lu,*) 'CHARAC: WITH STOCHASTIC DIFFUSION'
          WRITE(lu,*) '        AN ARGUMENT VISC MUST BE GIVEN'
          CALL plante(1)
          stop
        ELSE
!         HERE A DUMMY TARGET, WILL NOT BE USED
!         MAYBE NULLIFY WOULD BE BETTER ?
          avisc => ifamas
        ENDIF
      ENDIF
!
!-----------------------------------------------------------------------
!  TABLEAUX DE TRAVAIL PRIS DANS LE BLOC TB
!-----------------------------------------------------------------------
!
      IF(tb%N.GE.10) THEN
        t8 =>tb%ADR( 8)%P
        t9 =>tb%ADR( 9)%P
        t10=>tb%ADR(10)%P
        IF(optcha.GT.1) THEN
          t1 =>t1weak
          t2 =>t2weak
          t3 =>t3weak
          t4 =>t4weak
          t5 =>t5weak
          t6 =>t6weak
          t7 =>t7weak
          pt_shp=>shpwea%R
          pt_shz=>shzwea%R
          pt_shpbuf=>shpbuf
          pt_ftild=>ftild_weak
        ELSE
          t1 =>tb%ADR(1)%P
          t2 =>tb%ADR(2)%P
          t3 =>tb%ADR(3)%P
          t4 =>tb%ADR(4)%P
          t5 =>tb%ADR(5)%P
          t6 =>tb%ADR(6)%P
          t7 =>tb%ADR(7)%P
          pt_shp=>shp%R
          pt_shz=>shz%R
          pt_shpbuf=>mesh%M%X
          pt_ftild=>ftild
        ENDIF
      ELSE
        WRITE(lu,*) 'SIZE OF BLOCK TB:',tb%N
        WRITE(lu,*) 'TOO SMALL IN CHARAC'
        WRITE(lu,*) '10 REQUESTED'
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!  DEPLOIEMENT DE LA STRUCTURE DE MAILLAGE
!-----------------------------------------------------------------------
!
      npoin = mesh%NPOIN
      ielmu = uconv%ELM
!
!     PREPARING WORK ARRAYS
!
!     THE OFF-DIAGONAL TERMS OF WORK MATRIX IN MESH WILL BE USED AS
!     SHPBUF(3,SIZEBUF)
!
      IF(optcha.GT.1) THEN
        sizebuf=ng
      ELSE
        sizebuf=(mesh%M%X%MAXDIM1*mesh%M%X%MAXDIM2)/3
!       T7 WILL BE USED AS SHZBUF(SIZEBUF)
        sizebuf=min(sizebuf,t7%MAXDIM1)
!       IT3 WILL BE USED AS ELTBUF
        sizebuf=min(sizebuf,npoin)
      ENDIF
!
!-----------------------------------------------------------------------
!  ARE THERE QUADRATIC OR QUASI-BUBBLE VARIABLES ?
!  AND COMPUTATION OF LARGEST NUMBER OF POINTS
!-----------------------------------------------------------------------
!
      quad=.false.
      quab=.false.
      npt=0
      IF(fn%TYPE.EQ.4) THEN
        DO i=1,fn%N
          IF(fn%ADR(i)%P%ELM.EQ.12) quab = .true.
          IF(fn%ADR(i)%P%ELM.EQ.13) quad = .true.
          npt=max(npt,fn%ADR(i)%P%DIM1)
        ENDDO
      ELSEIF(fn%TYPE.EQ.2) THEN
        IF(fn%ELM.EQ.12) quab = .true.
        IF(fn%ELM.EQ.13) quad = .true.
        npt=max(npt,fn%DIM1)
      ENDIF
      IF(quab.AND.quad) THEN
        WRITE(lu,*) 'CHARAC: QUADRATIC AND QUASI-BUBBLE CANNOT BE MIXED'
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!     CHECKING SHP SIZE (ONCE A BUG...)
!-----------------------------------------------------------------------
!
      IF(3*npt.GT.shp%MAXDIM1*shp%MAXDIM2) THEN
        WRITE(lu,*) 'SIZE OF SHP:',shp%MAXDIM1*shp%MAXDIM2
        WRITE(lu,*) 'TOO SMALL IN CHARAC, ',3*npt
        WRITE(lu,*) 'REQUESTED'
        CALL plante(1)
        stop
      ENDIF
!
      IF(max(npoin2,npt).GT.t2%MAXDIM1) THEN
        WRITE(lu,*) 'SIZE OF T2:',t2%MAXDIM1
        WRITE(lu,*) 'TOO SMALL IN CHARAC, ',max(npoin2,npt)
        WRITE(lu,*) 'REQUESTED'
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!  APPEL DE SCARACT
!-----------------------------------------------------------------------
!
      IF(msk) THEN
!       APPEL AVEC IFAMAS
        ifa=>ifamas%I
      ELSE
!       APPEL AVEC IFABOR
        ifa=>mesh%IFABOR%I
      ENDIF
!
!     STARTING X AND Y OF POINTS (T1=XCONV AND T2=YCONV)
!
      IF(optcha.GT.1) THEN
!
        CALL char_gauss(t1%R,t2%R,t3%R,shpwea%R,shzwea%R,elt,eta,
     &                  mesh%X%R,mesh%Y%R,
     &                  ikle2%I,npoin2,nelem2,nelmax2,ng,ngauss,ielm,
     &                  nplan,zstar%R)
        IF(ielm.EQ.11) THEN
          dim1f=ng
        ELSEIF(ielm.EQ.41) THEN
          dim1f=npoin2
        ENDIF
        npt=ng
!
      ELSE
!
        CALL os('X=Y     ',x=t1,y=mesh%X)
        CALL os('X=Y     ',x=t2,y=mesh%Y)
!
!       IELM MUST BE INTENT(INOUT) BECAUSE IT IS SUCH IN CHGDIS
        IF(quad) THEN
          CALL chgdis(t1,ielm,13,mesh)
          CALL chgdis(t2,ielm,13,mesh)
        ELSEIF(quab) THEN
          CALL chgdis(t1,ielm,12,mesh)
          CALL chgdis(t2,ielm,12,mesh)
        ENDIF
!
        IF(ielm.EQ.11) THEN
          CALL gtsh11(shp%R,elt,ikle2%I,mesh%ELTCAR%I,npoin2,
     &                nelem2,nelmax2,mesh%NSEG,quab,quad)
          dim1f=npt
        ELSEIF(ielm.EQ.41) THEN
!         STARTING Z OF POINTS (T3=ZCONV)
          DO iplan=1,nplan
            DO i= (iplan-1)*npoin2+1,iplan*npoin2
              t3%R(i)=zstar%R(iplan)
            ENDDO
          ENDDO
!         IN 4D, STARTING F OF POINTS (T9=FCONV)
          IF(aya4d) THEN
            CALL ov('X=C     ', x=t9%R, c=freq%R(jf),
     &              dim1=npoin2*nplan)
          ENDIF
          CALL gtsh41(shp%R,shz%R,shf%R,wconv%R,frconv%R,
     &                elt,eta,fre,ikle2%I,mesh%ELTCAR%I,
     &                npoin2,nelmax2,nplan,jf,nf,aya4d)
          dim1f=npoin2
        ELSE
          WRITE(lu,*) 'ELEMENT NOT IMPLEMENTED IN CHARAC: ',ielm
          CALL plante(1)
          stop
        ENDIF
!
      ENDIF
!
      CALL scaract(fn,pt_ftild,uconv%R,vconv%R,wconv%R,frconv%R,
     &             mesh%X%R,mesh%Y%R,zstar%R,freq%R,
!                  XCONV YCONV ZCONV FCONV DX   DY   DZ   DF
     &             t1%R ,t2%R ,t3%R ,t9%R ,t4%R,t5%R,t6%R,t8%R,
!                           SHP
     &             mesh%Z%R,pt_shp,pt_shz,shf%R,
     &             surdet2%R,dt,ikle2%I,ifa,elt,
     &             eta,fre,it3,isub,
     &             ielm,ielmu,nelem2,nelmax2,nomb,npoin2,
     &             3,nrk,nplan,nf,mesh,npt,dim1f,-1,
!                  SHPBUF      SHZBUF      SHFBUF
     &             pt_shpbuf%R,t7%R       ,t10%R,frebuf,sizebuf,
     &             apost,aperio,aya4d,asigma,astocha,avisc)
!
      IF(optcha.GT.1) THEN
!
        IF(nomb.GT.0) THEN
          IF(ftild%TYPE.EQ.2) THEN
            CALL char_weak(ftild,ftild_weak,mesh%SURFAC%R,ikle2%I,
     &                     npoin2,nelem2,nelmax2,ng,ngauss,
     &                     mesh,t2,tb,agglo,ielm,nplan,mesh%Z%R,
     &                     rhs,am1,slv,unsv,listin,.true.)
          ELSEIF(ftild%TYPE.EQ.4) THEN
            DO i=1,nomb
              CALL char_weak(ftild%ADR(i)%P,ftild_weak%ADR(i)%P,
     &                       mesh%SURFAC%R,ikle2%I,
     &                       npoin2,nelem2,nelmax2,ng,ngauss,
     &                       mesh,t2,tb,agglo,ielm,nplan,mesh%Z%R,
     &                       rhs,am1,slv,unsv,listin,.true.)
            ENDDO
          ENDIF
        ENDIF
!
      ELSE
!
!     PARALLEL COMMUNICATION
!
        IF(ncsize.GT.1.AND.nomb.GT.0) THEN
          IF(ftild%TYPE.EQ.2) THEN
            CALL parcom(ftild,1,mesh)
          ELSEIF(ftild%TYPE.EQ.4) THEN
            DO i=1,nomb
              CALL parcom(ftild%ADR(i)%P,1,mesh)
            ENDDO
          ENDIF
        ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END