C++
C++ DENSITY CALCULATION
C++
      SUBROUTINE DENSIT(NO1,NO2,VC,CUX,CUY,CUZ,G1X)

      IMPLICIT NONE

!     Les parametres du programme
      INCLUDE  "PARAM"
      PARAMETER (FIELD1=FIELD3+2)
    
!     Les arguments transmis
      integer no1, no2
      double precision vc
      double precision CUX(FIELD1,FIELD2), CUY(FIELD1,FIELD2)
      double precision CUZ(FIELD1,FIELD2), G1X(FIELD1,FIELD2)

!     Les communs utilises
      double precision x, y
      COMMON/PARTI/X(PART),Y(PART)
      double precision px, py, pz, psa
      COMMON/VITE/PX(PART),PY(PART),PZ(PART),PSA(PART)
      integer no,ncx,ncy,lvgc
      COMMON/CONST/NO,NCX,NCY,LVGC
      double precision rmass, rte
      COMMON/CONST1/RMASS,RTE
      double precision bncxr,bncxl,bncxl2,bncxr2,
     & ancyl,ancyr,boundy,boundx
      COMMON/BOUND/BNCXR,BNCXL,BNCXL2,BNCXR2,ANCYL,ANCYR,BOUNDY,BOUNDX
      integer ir, il, jr, jl
      COMMON/PERIOD/IR(FIELD3),IL(FIELD3),JR(FIELD2),JL(FIELD2)

!     Tableaux sur les grilles
      integer ngrid
      parameter (ngrid=35)
      real cux_g(FIELD1,FIELD2,ngrid),cuy_g(FIELD1,FIELD2,ngrid),
     &     cuz_g(FIELD1,FIELD2,ngrid),g1x_g(FIELD1,FIELD2,ngrid)
      integer ig

!     Variables locales
      double precision vx, vy, vz
      double precision vax, vay
      double precision dx, dy, dxx, dyy, dxy, dyx, dzx, dzy
      integer lx, ly, lxl, lxr, lyl, lyr

!     Indices de boucles
      integer i1, i2, i
      integer m1, m2

!     Mesure du temps
      external second

C++
C++  FOR SPECULARLY REFLECTION
C++
      logical resist
      COMMON/RESI/RESIST

!     Mise a zero initiale
      cux_g = 0
      cuy_g = 0
      cuz_g = 0
      g1x_g = 0

      m1 = (no1-1)*2048+1
      m2 = no2*2048

      do i = m1,m2
C++
C++   REFLECT PARTICLES AT WALL
C++
          VAX=VC*PX(I)/PSA(I)
          VAY=VC*PY(I)/PSA(I)
          X(I)=X(I)+VAX
          Y(I)=Y(I)+VAY
          IF(Y(I).LT.ANCYL) Y(I)=Y(I)+BOUNDY
          IF(Y(I).GE.ANCYR) Y(I)=Y(I)-BOUNDY
C          IF(X(I).LT.ANCXL) X(I)=X(I)+BOUNDX
C          IF(X(I).GE.ANCXR) X(I)=X(I)-BOUNDX

           IF (X(I).LE.BNCXL) THEN
             X(I)=BNCXL2-X(I)
             IF (RESIST) THEN
               PX(I)=-PX(I)
               PY(I)=-PY(I)
             ELSE
               PX(I)=-PX(I)
               PZ(I)=-PZ(I)
              ENDIF
           ELSE IF(X(I).GE.BNCXR) THEN
             X(I)=BNCXR2-X(I)
              IF (RESIST) THEN
                PX(I)=-PX(I)
                PY(I)=-PY(I)
              ELSE
                PX(I)=-PX(I)
                PZ(I)=-PZ(I)
              ENDIF
            ENDIF
C++
C++ END TEST
C++
        end do

      i1 = m1
      DO WHILE (i1.le.m2)
        i2 = min(i1+ngrid-1,m2)

CDIR$ IVDEP
CDIR$ SHORTLOOP
        DO I=I1,I2

          IG = I - I1 + 1

          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)

          VX = PX(I)/PSA(I)
          VY = PY(I)/PSA(I)
          VZ = PZ(I)/PSA(I)

          DXX=VX*DX
          DXY=VX*DY
          DYY=VY*DY
          DYX=VY*DX
          DZX=VZ*DX
          DZY=VZ*DY

          CUX_G(LX,LY,IG)  = CUX_G(LX,LY,IG)  + VX
          CUX_G(LX,LYR,IG) = CUX_G(LX,LYR,IG) + DXY
          CUX_G(LX,LYL,IG) = CUX_G(LX,LYL,IG) - DXY
          CUX_G(LXL,LY,IG) = CUX_G(LXL,LY,IG) - DXX
          CUX_G(LXR,LY,IG) = CUX_G(LXR,LY,IG) + DXX
          CUY_G(LX,LY,IG)  = CUY_G(LX,LY,IG)  + VY
          CUY_G(LX,LYR,IG) = CUY_G(LX,LYR,IG) + DYY
          CUY_G(LX,LYL,IG) = CUY_G(LX,LYL,IG) - DYY
          CUY_G(LXL,LY,IG) = CUY_G(LXL,LY,IG) - DYX
          CUY_G(LXR,LY,IG) = CUY_G(LXR,LY,IG) + DYX
          CUZ_G(LX,LY,IG)  = CUZ_G(LX,LY,IG)  + VZ
          CUZ_G(LX,LYR,IG) = CUZ_G(LX,LYR,IG) + DZY
          CUZ_G(LX,LYL,IG) = CUZ_G(LX,LYL,IG) - DZY
          CUZ_G(LXL,LY,IG) = CUZ_G(LXL,LY,IG) - DZX
          CUZ_G(LXR,LY,IG) = CUZ_G(LXR,LY,IG) + DZX

        END DO

        I1 = I2 + 1

      END DO

      do i = m1,m2
C++
C++   REFLECT PARTICLES AT WALL
C++
          VAX=VC*PX(I)/PSA(I)
          VAY=VC*PY(I)/PSA(I)
          X(I)=X(I)+VAX
          Y(I)=Y(I)+VAY
          IF(Y(I).LT.ANCYL) Y(I)=Y(I)+BOUNDY
          IF(Y(I).GE.ANCYR) Y(I)=Y(I)-BOUNDY
C          IF(X(I).LT.ANCXL) X(I)=X(I)+BOUNDX
C          IF(X(I).GE.ANCXR) X(I)=X(I)-BOUNDX

           IF (X(I).LE.BNCXL) THEN
             X(I)=BNCXL2-X(I)
             IF (RESIST) THEN
               PX(I)=-PX(I)
               PY(I)=-PY(I)
             ELSE
               PX(I)=-PX(I)
               PZ(I)=-PZ(I)
             ENDIF
           ELSE IF(X(I).GE.BNCXR) THEN
             X(I)=BNCXR2-X(I)
             IF (RESIST) THEN
               PX(I)=-PX(I)
               PY(I)=-PY(I)
             ELSE
               PX(I)=-PX(I)
               PZ(I)=-PZ(I)
             ENDIF
           ENDIF
C++
C++ END TEST
C++
      end do

      i1 = m1
      DO WHILE (i1.le.m2)
        i2 = min(i1+ngrid-1,m2)

CDIR$ IVDEP
CDIR$ SHORTLOOP
        DO I=I1,I2

          IG = I - I1 + 1

          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)

          G1X_G(LX,LY,IG)  = G1X_G(LX,LY,IG)  + 1.
          G1X_G(LX,LYR,IG) = G1X_G(LX,LYR,IG) + DY
          G1X_G(LX,LYL,IG) = G1X_G(LX,LYL,IG) - DY
          G1X_G(LXL,LY,IG) = G1X_G(LXL,LY,IG) - DX
          G1X_G(LXR,LY,IG) = G1X_G(LXR,LY,IG) + DX

        END DO

        I1 = I2 + 1

      END DO

!     On recolle les morceaux
      do ig=1,ngrid
        do i2=1,field2
          do i1=1,field1
            cux(i1,i2) = cux(i1,i2) + cux_g(i1,i2,ig)
            cuy(i1,i2) = cuy(i1,i2) + cuy_g(i1,i2,ig)
            cuz(i1,i2) = cuz(i1,i2) + cuz_g(i1,i2,ig)
            g1x(i1,i2) = g1x(i1,i2) + g1x_g(i1,i2,ig)
          end do
        end do
      end do

      END