*VERSION D<FINITIVE : - M<THODE D'ADAM
*                     - NOTION DE GROUPES
*                     - TRAVAIL SUR ZONE UTILE UNIQUEMEMT
*                     - I2-I1+1 PEUT NE PAS ETRE MULTIPLE DE NGRP
*                     - TRANSMISSION PAR ARGUMENTS
C *****************************
C  DIAGNOSTICS
C *****************************
C++
C++ DENSITY CALCULATION
C++
      SUBROUTINE DENSIT(I1,I2,VC,CUX,CUY,CUZ,G1X)

C--      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

        INCLUDE  "PARAM"
    
      PARAMETER (FIELD1=FIELD3+2)
      COMMON/PARTI/X(PART),Y(PART)
      COMMON/VITE/PX(PART),PY(PART),PZ(PART),PSA(PART)
      COMMON/CONST/NO,NCX,NCY,LVGC
      COMMON/CONST1/RMASS,RTE
      COMMON/BOUND/BNCXR,BNCXL,BNCXL2,BNCXR2,ANCYL,ANCYR,BOUNDY,BOUNDX
      COMMON/PERIOD/IR(FIELD3),IL(FIELD3),JR(FIELD2),JL(FIELD2)
*
      DIMENSION CUX(FIELD1,FIELD2),CUY(FIELD1,FIELD2),CUZ(FIELD1,FIELD2)
      DIMENSION G1X(FIELD1,FIELD2)
C++
C++  FOR SPECULARLY REFLECTION
C++
      COMMON/RESI/RESIST
      LOGICAL RESIST
*
*     NGRP DONNE LE NOMBRE DE GROUPES UTILIS{
*     ON TRAVAILLERA SUR NGRP GRILLES IND{PENDANTES QUE L'ON SOMMERA
*     @ LA FIN
*     IDEB ET IFIN MARQUENT LE D{BUT DE LA ZONE UTILE DE CUX ... G1X
*     CDCOPY D{FINIT LE TABLEAU DE TRAVAIL

      PARAMETER (IDEB=LPLXBD-1)

CCC parametres VP200      PARAMETER(NGRP=12,IFIN=FIELD3,LDCDC=IFIN-IDEB+1)
      PARAMETER(NGRP=32,IFIN=FIELD3,LDCDC=IFIN-IDEB+1)
      DOUBLE PRECISION CDCOPY(IDEB:IFIN,FIELD2,4,NGRP)
      DOUBLE PRECISION CDCEQU(IDEB:IDEB+LDCDC*FIELD2*4*NGRP)
      EQUIVALENCE (CDCOPY(IDEB,1,1,1),CDCEQU(IDEB))
*
*     LONVEC EST LE PARAM}TRE DE LA M{THODE D'ADAM
*     (PUISQU'ON VA MODIFIER CHACUN DES 4 TABLEAUX CUX ... G1X
*     EN 5 POINTS (I,J) , (I-1,J) ETC .
*     INDICE STOCKERA PROVISOIREMENT L'INDICE @ MODIFIER ET AJOUT
*     STOCKERA LE NOMBRE @ AJOUTER AU TERME AYANT CET INDICE
*     (M{THODE D'ADAM)
*     LA M{THODE D'ADAM MODIFI{E CONSISTE @ VECTORISER @ LA FOIS
*     SUR NGRP ET SUR LONVEC, SOIT SUR UNE LONGUEUR DE VECTEUR
*     DE NGRP*LONVEC. D'O: LES {QUIVALENCES UTILIS{ES
CCCCC      PARAMETER(LREGIS=2048,LONVEC=20)
      PARAMETER(LREGIS=2048,LONVEC=20)
*     SURDIMENSIONNEMENT 1ER INDICE POUR RAISONS DE STRIDE
      PARAMETER(LDIND=LREGIS+2,LDAJ=LREGIS+1)
      DIMENSION INDICE(LDIND,LONVEC,NGRP),AJOUT(LDAJ,LONVEC,NGRP)
      DIMENSION INDEQU(LDIND,LONVEC*NGRP),AJEQU(LDAJ,LONVEC*NGRP)
      EQUIVALENCE (INDICE(1,1,1),INDEQU(1,1))
      EQUIVALENCE (AJOUT (1,1,1),AJEQU (1,1))
*
*     TERME D'INDICE DU AU TABLEAU UTILIS{ DANS CD (CUX, ... G1X)
      INDCUX = 0
      INDCUY =   LDCDC*FIELD2
      INDCUZ = 2*LDCDC*FIELD2
      INDG1X = 3*LDCDC*FIELD2
*
C++++++
C 
C       LXMAX=0
C       LXMIN=9000
*     INITIALISATION ZONE DE TRAVAIL
      DO 7 K = IDEB,IDEB+LDCDC*FIELD2*4*NGRP
        CDCEQU(K) = 0.D0
  7   CONTINUE
*
*     AU CAS O: I2-I1+1 NE SERAIT PAS MULTIPLE DE LREGIS * NGRP
      IQUOT = (I2-I1)/NGRP
      IREST = MOD(I2-I1,NGRP)
*
      DO 17 JQUOT = 1,IQUOT+1
*
        M = I1 + (JQUOT-1) * NGRP
*
        IF (JQUOT.NE.(IQUOT+1)) THEN
          JREST = NGRP
        ELSE
          JREST = IREST+1
*         PR{CAUTIONS INITIALISATION PARTIE INUTILIS{E D'AJOUT
          IF (JREST.LT.NGRP) THEN
            DO 22 NG=JREST+1,NGRP
              DO 21 K = 1,LONVEC
                DO 20 J = 1,LREGIS
                  AJOUT(J,K,NG) = 0.D0
   20           CONTINUE
   21         CONTINUE
   22       CONTINUE
          ENDIF
        ENDIF
*
        DO 18 NG=1,JREST
*
*       TERME D'INDICE DU AU GROUPE
        INDGRP = (NG-1)*LDCDC*FIELD2*4
*
        DO 19 J=1,LREGIS
*
*         INDICE DE LA PARTICULE UTILIS{E
          I=(M-1)*LREGIS+(NG-1)*LREGIS+J
*
          VX=PX(I)/PSA(I)
          VY=PY(I)/PSA(I)
          VZ=PZ(I)/PSA(I)
          VAX=VC*VX
          VAY=VC*VY
          X(I)=X(I)+VAX
          Y(I)=Y(I)+VAY
C--       IF(X(I).LT.ANCXL) X(I)=X(I)+BOUNDX
C--       IF(X(I).GE.ANCXR) X(I)=X(I)-BOUNDX
CCCCC*VOCL STMT,IF(2)
CCC*VDIR IF(2)
          IF(Y(I).LT.ANCYL) Y(I)=Y(I)+BOUNDY
CCCCCCCCC*VOCL STMT,IF(2)
CCC*VDIR IF(2)
          IF(Y(I).GE.ANCYR) Y(I)=Y(I)-BOUNDY
C++
C++   REFLECT PARTICLES AT WALL
C++
CCCCCCCC*VOCL STMT,IF(2)
CCC*VDIR IF(2)
          IF (X(I).LE.BNCXL) THEN
            X(I)=BNCXL2-X(I)
            IF (RESIST) THEN
            PX(I)=-PX(I)
            PY(I)=-PY(I)
            VX=-VX
            VY=-VY
            ELSE
            PX(I)=-PX(I)
            PZ(I)=-PZ(I)
            VX=-VX
            VZ=-VZ
            ENDIF
          ELSE
CCCCCCCCCCC*VOCL STMT,IF(2)
CCCC*VDIR IF(2)
            IF(X(I).GE.BNCXR) THEN
              X(I)=BNCXR2-X(I)
              IF (RESIST) THEN
              PX(I)=-PX(I)
              PY(I)=-PY(I)
              VX=-VX
              VY=-VY
              ELSE
              PX(I)=-PX(I)
              PZ(I)=-PZ(I)
              VX=-VX
              VZ=-VZ
              ENDIF
            ENDIF
          ENDIF
C++
C++ END TEST
C++
          LX=X(I)+1.5
          DX=(X(I)-LX+1)*0.5
          LXL=IL(LX)
          LXR=IR(LX)
          LY=Y(I)+1.5
          DY=(Y(I)-LY+1)*0.5
          LYL=JL(LY)
          LYR=JR(LY)
*
          DXX=VX*DX
          DXY=VX*DY
          DYY=VY*DY
          DYX=VY*DX
          DZX=VZ*DX
          DZY=VZ*DY
*
*         TERMES D'INDICE DUS AUX POINTS DE GRILLE CONSID{R{S
          LXLY  = LX  + LDCDC*(LY -1)
          LXLYR = LX  + LDCDC*(LYR-1)
          LXLYL = LX  + LDCDC*(LYL-1)
          LXLLY = LXL + LDCDC*(LY -1)
          LXRLY = LXR + LDCDC*(LY -1)
*JM(
*         INDICE COMPLET SOMME DE TROIS TERMES
*         - LE TERME DU AUX POINTS DE GRILLE
*         - LE TERME DU AU TABLEAU CONSID{R{ (CUX ... G1X)
*         - LE TERME DU AU GROUPE CHOISI
          INDICE(J, 1,NG) =   LXLY  + INDCUX + INDGRP
          AJOUT (J, 1,NG) = VX
*JM-      CUX(LX,LY)=CUX(LX,LY)+VX
          INDICE(J, 2,NG) =   LXLYR + INDCUX + INDGRP
          AJOUT (J, 2,NG) = DXY
*JM-      CUX(LX,LYR)=CUX(LX,LYR)+DXY
          INDICE(J, 3,NG) =   LXLYL + INDCUX + INDGRP
          AJOUT (J, 3,NG) = - DXY
*JM-      CUX(LX,LYL)=CUX(LX,LYL)-DXY
          INDICE(J, 4,NG) =   LXLLY + INDCUX + INDGRP
          AJOUT (J, 4,NG) = - DXX
*JM-      CUX(LXL,LY)=CUX(LXL,LY)-DXX
          INDICE(J, 5,NG) =   LXRLY + INDCUX + INDGRP
          AJOUT (J, 5,NG) = DXX
*JM-      CUX(LXR,LY)=CUX(LXR,LY)+DXX
          INDICE(J, 6,NG) =   LXLY  + INDCUY + INDGRP
          AJOUT (J, 6,NG) = VY
*JM-      CUY(LX,LY)=CUY(LX,LY)+VY
          INDICE(J, 7,NG) =   LXLYR + INDCUY + INDGRP
          AJOUT (J, 7,NG) = DYY
*JM-      CUY(LX,LYR)=CUY(LX,LYR)+DYY
          INDICE(J, 8,NG) =   LXLYL + INDCUY + INDGRP
          AJOUT (J, 8,NG) = - DYY
*JM-      CUY(LX,LYL)=CUY(LX,LYL)-DYY
          INDICE(J, 9,NG) =   LXLLY + INDCUY + INDGRP
          AJOUT (J, 9,NG) = - DYX
*JM-      CUY(LXL,LY)=CUY(LXL,LY)-DYX
          INDICE(J,10,NG) =   LXRLY + INDCUY + INDGRP
          AJOUT (J,10,NG) = DYX
*JM-      CUY(LXR,LY)=CUY(LXR,LY)+DYX
          INDICE(J,11,NG) =   LXLY  + INDCUZ + INDGRP
          AJOUT (J,11,NG) = VZ
*JM-      CUZ(LX,LY)=CUZ(LX,LY)+VZ
          INDICE(J,12,NG) =   LXLYR + INDCUZ + INDGRP
          AJOUT (J,12,NG) = DZY
*JM-      CUZ(LX,LYR)=CUZ(LX,LYR)+DZY
          INDICE(J,13,NG) =   LXLYL + INDCUZ + INDGRP
          AJOUT (J,13,NG) = - DZY
*JM-      CUZ(LX,LYL)=CUZ(LX,LYL)-DZY
          INDICE(J,14,NG) =   LXLLY + INDCUZ + INDGRP
          AJOUT (J,14,NG) = - DZX
*JM-      CUZ(LXL,LY)=CUZ(LXL,LY)-DZX
          INDICE(J,15,NG) =   LXRLY + INDCUZ + INDGRP
          AJOUT (J,15,NG) = DZX
*JM-      CUZ(LXR,LY)=CUZ(LXR,LY)+DZX
*JM)
C++
C++  ADVANCE PARTICLES ONE HALF-STEP
C++  CALCULATE CURRENT
C++
          VAX=VC*VX
          VAY=VC*VY
          X(I)=X(I)+VAX
          Y(I)=Y(I)+VAY
C--       IF(X(I).LT.ANCXL) X(I)=X(I)+BOUNDX
C--       IF(X(I).GE.ANCXR) X(I)=X(I)-BOUNDX
CCCCCCCCCCCCCC*VOCL STMT,IF(2)
CCCCCCC*VDIR IF(2)
          IF(Y(I).LT.ANCYL) Y(I)=Y(I)+BOUNDY
CCCCCCCCCCCCCCC*VOCL STMT,IF(2)
CCCCCCCCCC*VDIR IF(2)
          IF(Y(I).GE.ANCYR) Y(I)=Y(I)-BOUNDY
C++
C++   REFLECT PARTICLES AT WALL
C++
CCCCCCCCCCCC*VOCL STMT,IF(2)
CCCCCCCCCC*VDIR IF(2)
          IF (X(I).LE.BNCXL) THEN
            X(I)=BNCXL2-X(I)
            IF (RESIST) THEN
            PX(I)=-PX(I)
            PY(I)=-PY(I)
            VX=-VX
            VY=-VY
            ELSE
            PX(I)=-PX(I)
            PZ(I)=-PZ(I)
            VX=-VX
            VZ=-VZ
            ENDIF
          ELSE
CCCCCCCCCC*VOCL STMT,IF(2)
CCCCCCCCCCCC*VDIR IF(2)
            IF(X(I).GE.BNCXR) THEN
              X(I)=BNCXR2-X(I)
              IF (RESIST) THEN
              PX(I)=-PX(I)
              PY(I)=-PY(I)
              VX=-VX
              VY=-VY
              ELSE
              PX(I)=-PX(I)
              PZ(I)=-PZ(I)
              VX=-VX
              VZ=-VZ
              ENDIF
            ENDIF
          ENDIF
C++
C++ END TEST
C++
          LX=X(I)+1.5
          DX=(X(I)-LX+1)*0.5
          LXL=IL(LX)
          LXR=IR(LX)
          LY=Y(I)+1.5
          DY=(Y(I)-LY+1)*0.5
          LYL=JL(LY)
          LYR=JR(LY)

C
C         FOR NON PERIODIC CONDITIONS
C          DXL = DX
C          DXR = DX

C           IF (LX .EQ. LPLXBD-1) THEN
C           DXL = -DX
C            LXL=LX
C            ENDIF
C            IF (LX .GE. NCX-1)  THEN
C            DXR = -DX
C            LXR=LX
C            ENDIF
*
*         RECALCUL DES TERMES DU AUX POINTS DE GRILLE
          LXLY  = LX  + LDCDC*(LY -1)
          LXLYR = LX  + LDCDC*(LYR-1)
          LXLYL = LX  + LDCDC*(LYL-1)
          LXLLY = LXL + LDCDC*(LY -1)
          LXRLY = LXR + LDCDC*(LY -1)
*JM(
          INDICE(J,16,NG) =   LXLY  + INDG1X + INDGRP
          AJOUT (J,16,NG) = 1.
*JM-      G1X(LX,LY)=G1X(LX,LY)+1.
          INDICE(J,17,NG) =   LXLYR + INDG1X + INDGRP
          AJOUT (J,17,NG) = DY
*JM-      G1X(LX,LYR)=G1X(LX,LYR)+DY
          INDICE(J,18,NG) =   LXLYL + INDG1X + INDGRP
          AJOUT (J,18,NG) = - DY
*JM-      G1X(LX,LYL)=G1X(LX,LYL)-DY
          INDICE(J,19,NG) =   LXLLY + INDG1X + INDGRP
           AJOUT (J,19,NG) = - DX
*JM-      G1X(LXL,LY)=G1X(LXL,LY)-DX
          INDICE(J,20,NG) =   LXRLY + INDG1X + INDGRP
           AJOUT (J,20,NG) =   DX
*JM-      G1X(LXR,LY)=G1X(LXR,LY)+DX
C         XBF(I)=XBF(I)+DT*VX
*
   19   CONTINUE
*       (FIN DE LA BOUCLE SUR LE PAQUET DE LREGIS PARTICULES)
   18   CONTINUE
*       (FIN DE LA BOUCLE SUR LES GROUPES DE LREGIS PARTICULES)
*
*       ON APPORTE LES CONTRIBUTIONS CALCUL{ES AU TABLEAU DE TRAVAIL
*       ON NE PEUT PAS VECTORISER SUR J, MAIS ON PEUT VECTORISER SUR K
*       (ON A TOUT FAIT POUR)
        DO 24 J=1,LREGIS
CCCCCCCCCCCC*VOCL     LOOP,NOVREC
!CDIR     NODEP
CCCCCCCCCCCCCC*VDIR NODEP
          DO 23 K=1,LONVEC*NGRP
            CDCEQU(INDEQU(J,K)) = CDCEQU(INDEQU(J,K)) + AJEQU(J,K)
   23     CONTINUE
   24   CONTINUE
*
   17 CONTINUE
*     (FIN DE LA BOUCLE SUR LES PAQUETS DE NGRP*LREGIS PARTICULES)
*
*     ADDITION DES NGRP COLONNES DE CDCOPY POUR FORMER CUX, CUY ETC.
      DO 36 NG=1,NGRP
        DO 35 J=1,FIELD2
          DO 34 I=IDEB,IFIN
            CUX(I,J) = CUX(I,J) + CDCOPY(I,J,1,NG)
            CUY(I,J) = CUY(I,J) + CDCOPY(I,J,2,NG)
            CUZ(I,J) = CUZ(I,J) + CDCOPY(I,J,3,NG)
            G1X(I,J) = G1X(I,J) + CDCOPY(I,J,4,NG)
 34       CONTINUE
 35     CONTINUE
 36   CONTINUE
*
*JM)
      RETURN
      END
