Sources Esope | Cast3M

sautpa
C SAUTPA    SOURCE    CHAT      05/01/13    03:11:37     5004
      SUBROUTINE SAUTPA(EPSILO,ICHGZ,JREBO,XIREB,XNREB
     $     ,NDIM,ITY1,ITYG,NOEL1,NOUN1,DIAM,DIFL,DIFT,TDEP,XDEP2
     $     ,UELEM,DTREEL,UREEL,U,TARI,XCON2,XARI2,IZSH,Z
     $     ,IREBCO,IREBDI,NTEST)
*
*************************************************************************
*** SP 'SAUTPA' : fait avancer la particule apres un pas de tps 'DTREEL'
*** par convection (analytique ou explicite) + dispersion (explicite).
***
*** APPELES 1 = 'TDRAND'
*** APPELES 2 = 'REEREF', 'REFREE', 'JCBIEN', 'SAUCO1', SAUCO2',
***             'SAUDIF', 'NORMVI'
***
*** E = 'EPSILO' marge relative acceptée % 0 (cause machine)
***     'ICHGZ' vaut 1 si saut précédent effectif, 0 sinon
***     'JREBO' n° local face impermeable ou se trouve particule, -1 sinon
***     'XIREB' pt d'impact sur la face impermeable
***     'XNREB' vecteur normal à la face impermeable
***     'NDIM' dimension de l'espace
***     'ITY1' entier caracterisant le type de l'element
***     'ITYG' entier caracterisant la geometrie de l'element
***     'NOEL1' nombre de noeuds de l'element traversé
***     'NOUN1' nombre de flux de l'element traversé
***     'DIAM' "longueur caracteristique" de l'element traversé
***     'DIFL' coeff dispersion longitudinal ds element considere
***     'DIFT' coeff dispersion transversal ds element considere
***     'TDEP' tps reel depart avant l'avancee de la particule
***     'XDEP2' coord reelles de depart particule
***     'UELEM' valeurs des flux aux faces
***     'DTREEL' pas de tps réel pour le calcul avancée particule
***
*** S = 'UREEL' vecteur "vitesse convective" ds element reel
***     'U' vecteur "vitesse diffusive" ds element reel
***     'TARI' tps reel d'arrivee apres avancee particule
***     'XCON2' coord reelles arrivee particule apres convection
***     'XARI2' coord reelles arrivee particule apres diffusion
***     'NTEST' vaut 0 si Jacobien nul pour approximation vitesse, 0 sinon
***
*** E/S = 'IZSH' segment content pour l'elemt considere au pt depart :
***             'XYZL' coordonnees reelles des noeuds (E)
***             'SHP' fonctions de forme et derivees ds elemt ref (S)
***             'SHY' fonctions de base et derivees ds elemt ref (S)
***       'Z' vecteur aleatoire entre -1 et 1 pour le saut diffusif
***       'IREBCO' vaut 1 si rebond lors de la convection, 0 sinon
***       'IREBDI' vaut 1 si rebond lors de la diffusion, 0 sinon
***
***               Auteur Cyril Nou
*************************************************************************
*
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      SEGMENT IZSH
          REAL*8 SHP(6,MNO9),SHY(12,MNO9),XYZL(3,MNO9)
      ENDSEGMENT
      DIMENSION XDEP(3),XDEP2(3),XCON2(3),XARI2(3)
      DIMENSION UELEM(NOUN1),UREEL(3),U(3),XJAC(3,3),Z(3)
      DIMENSION XIREB(3),XNREB(3),XITEST(3),XNTEST(3),X(3)
 
C     write(6,*)' uelem ', (UELEM(ii),ii=1,NOUN1)
******************************
*** INITIALISATION DU SAUT ***
******************************
 
***   calcul position de départ de référence
      CALL REEREF(NDIM,ITY1,NOEL1,IZSH,XDEP2,XDEP)
***   initialisation à position de départ
      TARI=TDEP
      IREBCO=0
      IREBDI=0
      JFTEST=0
      DO 10 I=1,NDIM
         XCON2(I)=XDEP2(I)
         XARI2(I)=XDEP2(I)
 10   CONTINUE
 
**********************
*** SAUT CONVECTIF ***
**********************
 
***   cas ou le saut convectif est calculé explicitement
      CALL SAUCON(NDIM,ITY1,ITYG,NOEL1,NOUN1
     $     ,DIAM,UELEM,XDEP,XDEP2,DTREEL,XCON2,IZSH,UREEL,NTEST)
***   cas ou Jacobien nul lors approximation vitesse
      IF (NTEST.EQ.0) RETURN
***   calcul de la norme de la vitesse
      CALL NORMVI(NDIM,UREEL,XNORM,UXY)
********************************************************
*** CAS OU TRAVERSEE FACE IMPERMEABLE PDT CONVECTION ***
********************************************************
 
      CALL LIEUPT(0.D0,NDIM,ITYG,XCON2,IZSH,ITEST)
      IF (ITEST.EQ.0) THEN
         CALL CHOINT(NDIM,ITYG,XDEP2
     $        ,XDEP2,XCON2,IZSH,XITEST,XNTEST,KTEST,JFTEST)
         IF ((KTEST.EQ.2).AND.(JFTEST.EQ.JREBO)) THEN
            IREBCO=1
            COEFC=SCVECT(UREEL,XNREB,NDIM)
            COEFC=COEFC*DTREEL
            DO 30 I=1,NDIM
               XCON2(I)=XCON2(I)-IREBCO*COEFC*XNREB(I)
 30         CONTINUE
         ENDIF
      ENDIF
 
*********************************
*** SAUT DIFFUSIF (EXPLICITE) ***
*********************************
 
      IF (ICHGZ.EQ.1) THEN
***      on retire vecteur aleatoire car saut précédent effectif
         DO 20 I=1,3
            CALL TDRAND(Z(I))
            Z(I)=(2.D0*Z(I))-1.D0
 20      CONTINUE
      ENDIF
      CALL SAUDIF(NDIM,Z,XCON2,UREEL
     $     ,DTREEL,XNORM,UXY,DIFL,DIFT,U,XARI2)
 
*******************************************************
*** CAS OU TRAVERSEE FACE IMPERMEABLE PDT DIFFUSION ***
*******************************************************
 
      CALL LIEUPT(0.D0,NDIM,ITYG,XARI2,IZSH,JTEST)
      IF (JTEST.EQ.0) THEN
         CALL CHOINT(NDIM,ITYG,XCON2
     $        ,XCON2,XARI2,IZSH,XITEST,XNTEST,KTEST,JFTEST)
         IF ((KTEST.EQ.2).AND.(JFTEST.EQ.JREBO)) THEN
            IREBDI=1
            DO 40 I=1,NDIM
               XIREB(I)=XITEST(I)
               XNREB(I)=XNTEST(I)
 40         CONTINUE
            DO 50 I=1,NDIM
               X(I)=XARI2(I)-XIREB(I)
 50         CONTINUE
            COEFD=SCVECT(X,XNREB,NDIM)
            DO 60 I=1,NDIM
               XARI2(I)=XARI2(I)-2*IREBDI*COEFD*XNREB(I)
 60         CONTINUE
         ENDIF
      ENDIF
 
*******************************
*** TEMPS ARRIVEE PARTICULE ***
*******************************
 
      TARI=TDEP+DTREEL
      RETURN
      END