c@a
c@versb
C-----------------------------------------------------------------------
C
CVERS                  Code_Saturne version 1.3
C                      ------------------------
C
C     This file is part of the Code_Saturne Kernel, element of the
C     Code_Saturne CFD tool.
C
C     Copyright (C) 1998-2008 EDF S.A., France
C
C     contact: saturne-support@edf.fr
C
C     The Code_Saturne Kernel is free software; you can redistribute it
C     and/or modify it under the terms of the GNU General Public License
C     as published by the Free Software Foundation; either version 2 of
C     the License, or (at your option) any later version.
C
C     The Code_Saturne Kernel is distributed in the hope that it will be
C     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
C     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C     GNU General Public License for more details.
C
C     You should have received a copy of the GNU General Public License
C     along with the Code_Saturne Kernel; if not, write to the
C     Free Software Foundation, Inc.,
C     51 Franklin St, Fifth Floor,
C     Boston, MA  02110-1301  USA
C
C-----------------------------------------------------------------------
c@verse
                        SUBROUTINE CPPDF5
C                       *****************
C     -------------------------------------------------------------
     & ( NCELET , NCEL   ,
     &   F1M    , F2M    , F3M    , F4M    , F3P2M  , F4P2M  ,
     &   INDPDF ,
     &   SI7    , SI8    , SP2M   , F4I7   )
C     -------------------------------------------------------------
C***********************************************************************
C FONCTION :
C --------
c@foncb
CFONC
CFONC CALCUL DU TYPE DE PDF
CFONC   PDF CONJOINTE DEGENERRE EN UNE PDF 1D DE TYPE RECTANGLE - DIRAC
CFONC   SUPPORT (I7,I8) passant par M
CFONC
CFONC   --> RECONSTITUTION DE 5 MOMENTS
CFONC
c@fonce
C                             ARGUMENTS
c@argub
CARGU .______________.____._____.______________________________________.
CARGU !    NOM       !TYPE!MODE !                   ROLE               !
CARGU !______________!____!_____!______________________________________!
CARGU ! NCELET       ! E  !  -> ! NOMBRE D'ELEMENTS HALO COMPRIS       !
CARGU ! NCEL         ! E  !  -> ! NOMBRE D'ELEMENTS ACTIFS             !
CARGU ! F1M          ! TR !  -> ! MOYENNE DU TRACEUR 1 (MV CHZ +CO)    !
CARGU ! F2M          ! TR !  -> ! MOYENNE DU TRACEUR 2 (MV CXHY +CO)   !
CARGU ! F3M          ! TR !  -> ! MOYENNE DU TRACEUR 3 (CO C.HET)      !
CARGU ! F4M          ! TR !  -> ! MOYENNE DU TRACEUR 4 (AIR)           !
CARGU ! F3P2M        ! TR !  -> ! VARIANCE DU TRACEUR 3 (CO C.HET)     !
CARGU ! F4P2M        ! TR !  -> ! VARIANCE DU TRACEUR 4 (AIR)          !
CARGU ! INDPDF       ! TE !  <- ! PASSAGE PAR LES PDF                  !
CARGU ! SI7          ! TR !  <- ! ABSCISSE CURVILIGNE AU POINT I7      !
CARGU ! SI8          ! TR !  <- ! ABSCISSE CURVILIGNE AU POINT I8      !
CARGU ! SP2M         ! TR !  <- ! VARIANCE DROITE SUPPORT              !
CARGU ! F4I7         ! TR !  <- ! F4 AU POINT I7                       !
CARGU !______________!____!_____!______________________________________!
c@argue
C
c@commb
CCOMM                             COMMONS
CCOMM .______________.____._____.______________________________________.
CCOMM !    NOM       !TYPE!MODE !                   ROLE               !
CCOMM !______________!____!_____!______________________________________!
CCOMM !______________!____!_____!______________________________________!
c@comme
C
C     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
C            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
C     MODE : -> DONNEE, <- RESULTAT, <-> DONNEE MODIFIEE,
C            - TABLEAU DE TRAVAIL
C***********************************************************************
C
      IMPLICIT NONE
C
C***********************************************************************
C
      INCLUDE "paramx.h"
      INCLUDE "numvar.h"
      INCLUDE "optcal.h"
      INCLUDE "cstphy.h"
      INCLUDE "entsor.h"
      INCLUDE "cstnum.h"
      INCLUDE "ppppar.h"
      INCLUDE "ppthch.h"
      INCLUDE "coincl.h"
      INCLUDE "cpincl.h"
      INCLUDE "ppincl.h"
C
C***********************************************************************
C
C ARGUMENTS
C
      INTEGER          NCELET, NCEL
      INTEGER          INDPDF(NCELET)
C
      DOUBLE PRECISION F1M(NCELET), F2M(NCELET), F3M(NCELET)
      DOUBLE PRECISION F4M(NCELET), F3P2M(NCELET), F4P2M(NCELET)
      DOUBLE PRECISION SI7(NCELET), SI8(NCELET), SP2M(NCELET)
      DOUBLE PRECISION F4I7(NCELET)
C
C VARIABLES LOCALES
C
      INTEGER          IEL, II
      INTEGER          ITYP1, ITYPA, ITYP(2), ICPTS1
C
      DOUBLE PRECISION S2MAX, T1, T2
      DOUBLE PRECISION F4L5, RAV3S4
      DOUBLE PRECISION F1I7T(2), F2I7T(2), F3I7T(2), F4I7T(2), F4I8T(2)
      DOUBLE PRECISION SI8T(2), SI7T(2), SP2MT(2)
      DOUBLE PRECISION A, B, C, ANGALP, ANGGAM
      DOUBLE PRECISION F4I7TS, TETA(2)
      DOUBLE PRECISION XR1, F4I7T1, SI7T1, SI8T1, SP2MT1
      DOUBLE PRECISION XRA, F4I7TA, SI7TA, SI8TA, SP2MTA
      DOUBLE PRECISION F1I3, F2I3, F3I3, F1I5, F2I5, F3I5
      DOUBLE PRECISION XI4I5, XI4I3
C
C***********************************************************************
C
C=======================================================================
C 1. INITIALISATION
C=======================================================================
C
C --- Parametres relatifs a la variance T1 et a la moyenne T2
C
      T1 = 1.D-04
      T2 = 5.D-03
C
C --- Initialisation des tableaux
C
      DO IEL = 1, NCEL
        F4I7(IEL)   = 0.D0
        SI7(IEL)    = 0.D0
        SI8(IEL)    = 0.D0
        SP2M(IEL)   = 0.D0
        INDPDF(IEL) = 0.D0
      ENDDO
C
C
C=======================================================================
C 2. TYPE DE PDF
C    INDPDF = 0  : PAS DE PDF
C    INDPDF = 10 : SUPPORT (I7,I8) FORME AVEC TETA VARIABLE
C                  RECONSTITUTION DE 5 MOMENTS
C                    INDPDF = 11 : I7 APP. [L5,I3max]
C                    INDPDF = 12 : I7 APP. [I4,L3]
C                    INDPDF = 13 : I7 APP. [I4,L5]
C=======================================================================
C
      DO IEL = 1, NCEL
        IF ( F4P2M(IEL).GT.T1 .AND.
     &       F4M(IEL)  .GE.T2 .AND. F4M(IEL).LE.(1.D0-T2) ) THEN
          INDPDF(IEL) = 10
        ELSE
          INDPDF(IEL) = 0
        ENDIF
      ENDDO
C
C
C=======================================================================
C 3. CALCULS DES PARAMETRES DE LA PDF : SI7, SI8 et SP2M
C    CALCUL DE F4I7
C=======================================================================
C
!OCL SCALAR
C
      DO IEL = 1, NCEL
C
        IF ( INDPDF(IEL).EQ.10 ) THEN
C
          DO II = 1, 2
            F4I7T(II) = 0.D0
            SI7T(II)  = 0.D0
            SI8T(II)  = 0.D0
            SP2MT(II) = 0.D0
          ENDDO
C
C --> Calculs preliminaires
C
          RAV3S4 = F3P2M(IEL)/F4P2M(IEL)
C     Droite limite = (I5,M)
          F4L5 = F4M(IEL)/(F4M(IEL)+F3M(IEL))
C     Calcul de XI4I5
          F1I5 = F1M(IEL)/(F1M(IEL)+F2M(IEL))
          F2I5 = F2M(IEL)/(F1M(IEL)+F2M(IEL))
          F3I5 = 0.D0
          XI4I5 = SQRT (
     &            ( SQRT(6.D0)/2.D0*(0.D0-F1I5) +
     &              SQRT(6.D0)/4.D0*(0.D0-F2I5-F3I5)   )**2 +
     &            ( 3.D0*SQRT(2.D0)/4.D0*(0.D0-F2I5) +
     &                   SQRT(2.D0)/4.D0*(0.D0-F3I5)   )**2 +
     &            ( 0.D0-F3I5                          )**2 )
          F1I3 = 0.D0
          F2I3 = 0.D0
          F3I3 = 1.D0
C      Calcul de XI4I3
          XI4I3 = SQRT (
     &            ( SQRT(6.D0)/2.D0*(0.D0-F1I3) +
     &              SQRT(6.D0)/4.D0*(0.D0-F2I3-F3I3)  )**2 +
     &            ( 3.D0*SQRT(2.D0)/4.D0*(0.D0-F2I3) +
     &                   SQRT(2.D0)/4.D0*(0.D0-F3I3)  )**2 +
     &            ( 0.D0-F3I3                         )**2 )
          ANGALP = ACOS(XI4I3/(2.D0*XI4I5))
          ANGGAM = PI - 2.D0*ANGALP
C
C --> Determination des deux angles TETA
C
          A = 1.D0/(SIN(ANGALP)*XI4I3)
          B = 1.D0/XI4I5
          C = 1.D0/(XI4I5*TAN(ANGGAM))
          TETA(1) = ATAN(-SQRT(RAV3S4)*B/(A+SQRT(RAV3S4)*C))
          TETA(2) = ATAN( SQRT(RAV3S4)*B/(A-SQRT(RAV3S4)*C))
C
C --> Determination du type de PDF
C       - TYPE DE PDF : 13, I7 appartient au segment [I4,L5]
C       - TYPE DE PDF : 12, I7 appartient au segment [I4,L3]
C       - TYPE DE PDF : 11, I7 appartient au segment [L5,I3max]
C
          DO II = 1, 2
C
            A = SIN(TETA(II))/SIN(ANGALP)/XI4I3
            B = COS(TETA(II))/XI4I5+SIN(TETA(II))/XI4I5/TAN(ANGGAM)
            F4I7TS = (-A/B*F4M(IEL)+(1.D0-F3M(IEL)))/(1.D0-A/B)
C ---- Cette formulation permet de degenerer de facon naturelle
C      si RAV3S4 = 0.
C
            IF ( F4I7TS.LE.1.D0 .AND. F4I7TS.GT.F4L5 ) THEN
C
C ---- TYPE DE PDF : 13, I7 appartient au segment [I4,L5]
C                        I8 appartient au segment [I5,L4]
              ITYP(II) = 3
              F4I7T(II) = (-A/B*F4M(IEL)+(1.-F3M(IEL)))/(1.-A/B)
              F1I7T(II) = 0.D0
              F2I7T(II) = 0.D0
              F3I7T(II) = 1.D0 -F4I7T(II)
C
            ELSEIF ( F4I7TS.LE.F4L5 .AND.
     &               F4I7TS.GT.(1.D0-F3MAX) ) THEN
C
C ---- TYPE DE PDF : 11, I7 appartient au segment [L5,I3max]
C                        I8 appartient au segment [L3,I5]
              ITYP(II) = 1
              F4I7T(II) = (F4M(IEL)-(1.D0-F3M(IEL))*B/A)/(1.D0-B/A)
              F1I7T(II) = 0.D0
              F2I7T(II) = 0.D0
              F3I7T(II) = 1.D0 - F4I7T(II)
C
            ELSE
C
C TYPE DE PDF : 12, C appartient au segment [I4,L3]
C                   D appartient au segment [L4,I3max]
              ITYP(II) = 2
              F4I7T(II) = F3M(IEL)*B/A+F4M(IEL)
              F1I7T(II) = (1.D0-F4I7T(II))*F1M(IEL)/(F1M(IEL)+F2M(IEL))
              F2I7T(II) = (1.D0-F4I7T(II))*F2M(IEL)/(F1M(IEL)+F2M(IEL))
              F3I7T(II) = 0.D0
C
            ENDIF
C
            SI7T(II) = - SQRT (
     &  ( SQRT(6.D0)/2.D0*(F1M(IEL)-F1I7T(II)) +
     &    SQRT(6.D0)/4.D0*(F2M(IEL)+F3M(IEL)-F2I7T(II)-F3I7T(II)) )**2 +
     &  ( 3.D0*SQRT(2.D0)/4.D0*(F2M(IEL)-F2I7T(II)) +
     &         SQRT(2.D0)/4.D0*(F3M(IEL)-F3I7T(II))               )**2 +
     &  ( F3M(IEL)-F3I7T(II)                                      )**2 )
C
            F4I8T(II) =
     &  (1.D0-F3MAX)*(F4I7T(II)*(-F1M(IEL)-F2M(IEL)+1.D0)-F4M(IEL))
     &  / (F4I7T(II)-1.D0+F3M(IEL)+F3MAX*(F1M(IEL)+F2M(IEL)))
C
            SI8T(II) = SI7T(II)*
     &  (F4M(IEL)-F4I8T(II))/(F4M(IEL)-F4I7T(II))
C
            SP2MT(II) = F4P2M(IEL)/
     &  (F4M(IEL)-F4I7T(II))**2*(SI7T(II))**2
            S2MAX = -SI7T(II)*SI8T(II)
            IF (SP2MT(II).GT.S2MAX) SP2MT(II) = 0.D0
C
          ENDDO
C
C
C --> Vu que l'on a deux angles, on a deux supports possibles
C     On priviligie alors la forme 11 avec le support le plus long
C
          XRA     = ZERO
          SI7TA   = 0.D0
          SI8TA   = 0.D0
          SP2MTA  = 0.D0
          F4I7TA  = 0.D0
          ITYPA   = -10
C
          XR1     = ZERO
          SI7T1   = 0.D0
          SI8T1   = 0.D0
          SP2MT1  = 0.D0
          F4I7T1  = 0.D0
          ITYP1   = -10
          ICPTS1  = 0
C
          DO II = 1, 2
C ---- On priviligie la forme 11 avec le support le plus long
            IF ( (ITYP(II).EQ.1) ) THEN
              IF ( (SI8T(II)-SI7T(II)).GT.XR1 .AND.
     &             (SP2MT(II).GT.EPZERO)         ) THEN
                XR1     = (SI8T(II)-SI7T(II))
                F4I7T1  = F4I7T(II)
                ITYP1   = ITYP(II)
                SI7T1   = SI7T(II)
                SI8T1   = SI8T(II)
                SP2MT1  = SP2MT(II)
                ICPTS1  = ICPTS1 + 1
              ENDIF
            ELSE
              IF ( (SI8T(II)-SI7T(II)).GT.XRA .AND.
     &             (SP2MT(II).GT.EPZERO)         ) THEN
                XRA     = (SI8T(II)-SI7T(II))
                F4I7TA  = F4I7T(II)
                ITYPA   = ITYP(II)
                SI7TA   = SI7T(II)
                SI8TA   = SI8T(II)
                SP2MTA  = SP2MT(II)
              ENDIF
            ENDIF
          ENDDO
C
          IF ( ICPTS1.GE.1 ) THEN
            F4I7(IEL)   = F4I7T1
            SI7(IEL)    = SI7T1
            SI8(IEL)    = SI8T1
            SP2M(IEL)   = SP2MT1
            INDPDF(IEL) = INDPDF(IEL) + ITYP1
          ELSE
            F4I7(IEL)   = F4I7TA
            SI7(IEL)    = SI7TA
            SI8(IEL)    = SI8TA
            SP2M(IEL)   = SP2MTA
            INDPDF(IEL) = INDPDF(IEL) + ITYPA
          ENDIF
C
C
        ENDIF
C
      ENDDO
C
C
C
C----
C FIN
C----
C
      RETURN
      END
c@z
