C XFPA      SOURCE    CHAT      05/01/13    04:13:32     5004
      SUBROUTINE XFPA(ANU,YP,ROG,RAP,N,UET,VNORM,AK)
C*****************************************************************
C
C OBJET : calcul de la vitesse de depot d'aerosols AK dans chaque
C         element de la ligne paroi.
C
C SYNTAXE : AK = FPA ANU YP UET VNORM ROG RAP
C
C AK   : CHPOINT SCAL CENTRE (vitesse de depot)
C ANU  : FLOTTANT (viscosite)
C YP   : FLOTTANT (epaisseur de la couche limite)
C UET  : CHPOINT SCAL CENTRE (vitesse de frottement)
C VNORM : CHPOINT VECT FACE (normale a la paroi)
C ROG  : POINT (masse volumique * g)
C RAP  : FLOTTANT (rayon des particules)
C
C EST APPELEE PAR KFPA.SOR
C*****************************************************************
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION UET(*),VNORM(*),AK(*),ROG(*)

-INC PPARAM
-INC CCOPTIO
C
C
C
C
C ---------- parametres de la particule ------------------------------
C - coefficient de correction de Cunningham :
      CUN = 1.0D0 + 8.296D-8/RAP + 2.64D-8/(RAP*EXP(1.667D7*RAP))
C
C - masse volumique ROP :
      ROP = ROG(1)*ROG(1)
      DO 20 I=2,IDIM
            ROP = ROP + ROG(I)*ROG(I)
 20   CONTINUE
      ROP = SQRT(ROP)/9.81D0
C
C - vitesse de sedimentation incomplete VSI :
      VSI = 1.23D4*RAP*RAP*CUN
      VS  = VSI*ROP*9.81D0
C
C - temps de relaxation TAU :
      TAU = VSI*ROP
C
C - nombre de SCHMIDT :
      SCHMIDT = 5.62D21*ANU*ANU*RAP/CUN
C
C - constante de la vitesse de SAFFMAN (adimensionnelle) :
      CTSAF = 2.23D5*RAP*TAU
C
            DO 3453 IEL=1,N
C
C------- calcul de AKS = vitesse de sedimentation --------------------
C          ( d' apres l'orientation de l'element )
C
      PROSCA = ROG(1)*VNORM(IEL)
      DO 30 I=2,IDIM
            PROSCA = PROSCA + ROG(I)*VNORM(IEL+(I-1)*N)
 30   CONTINUE
      IF (PROSCA.LT.0.0D0) PROSCA = 0.0D0
      AKS = VSI*PROSCA
C
C
C------- calcul de AKDI = vitesse de depot par diffusion + inertie ----
C
      UIT = UET(IEL)
      RAPLUS = RAP*UIT/ANU
      TOPLUS = TAU*UIT*UIT/ANU
      YFIN = YP*UIT/ANU
C
C - initialisations et boucle :
      YPLUS1 = RAPLUS
      YPLUS2 = RAPLUS
      YPAS = 1.0D-4
      ANUT = 0.0D0
      ZETA1 = SCHMIDT
      ZETA2 = SCHMIDT
      ZETA = SCHMIDT*1.0D-4
      VPRIMF = 0.0D0
      VPRIMP1 = 0.0D0
      VPRIMP2 = 0.0D0
      TLPLUS = 0.0D0
      VMIGT = 0.0D0
      VSAFF = CTSAF
      VPLUS = 0.0D0
      CPLUS = 0.0D0
C
C
      DO 10   I = 1,200
              YPLUS2 = YPLUS2 + YPAS
         IF (YPLUS2.LE.YFIN) THEN
             IF (YPLUS2.LT.3) THEN
                 ANUT = (YPLUS2/1.115D1)**3
                 VSAFF = CTSAF/SQRT(1.0D0 + (YPLUS2/1.14D1)**2)
               ELSEIF (YPLUS2.LT.52) THEN
                 ANUT = (YPLUS2/1.14D1)**2 - 4.9774D-2
                 VSAFF = CTSAF/SQRT(1.0D0 + (YPLUS2/1.14D1)**2)
               ELSE
                 ANUT = 4.0D-1*YPLUS2
                 VSAFF = CTSAF/SQRT(1.0D0 + 4.0D-1*YPLUS2)
             ENDIF
C
             IF (YPLUS2.LT.30) THEN
                 VPRIM = 3.3D-2*(1.0D0 - EXP(-YPLUS2/3.837D0))
                 VPRIM = VPRIM*YPLUS2
                 VPR30 = VPRIM
               ELSE
                 IF (YFIN.EQ.3.0D1)   YFIN = 3.1D1
                 YPLUS3 = (YPLUS2-3.0D1)/(YFIN-3.0D1)
                 VPRIM = VPR30 - (VPR30 - 1.82D0)*YPLUS3
             ENDIF
                 VPRIMF = VPRIM*VPRIM
C
             ZETA2 = SCHMIDT/(1.0D0 + ANUT*SCHMIDT)
             ZETA = (ZETA1 + ZETA2)*YPAS/2.0D0
                   ZETA1 = ZETA2
C
             TLPLUS = ANUT/VPRIMF
             VPRIMP2 = VPRIMF/(1.0D0 + TOPLUS/TLPLUS)
             VMIGT = TOPLUS*(VPRIMP2 - VPRIMP1)/YPAS
                   VPRIMP1 = VPRIMP2
             VPLUS = VMIGT + VSAFF
             CPLUS = CPLUS + ZETA*(1.0D0 - CPLUS*VPLUS)
C
C       WRITE(6,*) 'CPLUS =' , CPLUS
C       WRITE(6,*) 'ZETA =' , ZETA
C       WRITE(6,*) 'VPLUS =' , VPLUS
C       WRITE(6,*) 'VPRIMP =' , VPRIMP2
C       WRITE(6,*) 'YPLUS =' , YPLUS2
             YPAS = (YPLUS2 - YPLUS1)*1.1D0
                   YPLUS1 = YPLUS2
         ENDIF
  10  CONTINUE
C
      AKPLUS = 1.0D0/CPLUS
      AKDI = AKPLUS*UIT
C
C
C------- calcul de AKT = vitesse de depot totale ----------------------
C
      AKT = AKS + AKDI
      AK(IEL) = AKT
C
C
 3453 CONTINUE
C
      RETURN
      END