!     ----------------------------------------------------------------
!     Compile as:  
!     f90 -O3 -n32 -r10000 -OPT:fast_io=ON mesh.f90
!     ----------------------------------------------------------------
      PROGRAM mesh
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MX=60,MY=30)
      INTEGER MVN(MX,MY)
      DOUBLE PRECISION MVX(MX,MY),MVY(MX,MY),MV2(MX,MY)
      DATA LP_DAT,LP_OUT,LP_TEMP,LP_DEN /25,26,27,28/      
      
      XMIN = 20.D0
      XMAX = 140.D0
      YMIN = 0.D0
      YMAX = 50.D0
      
      BINX = (XMAX-XMIN)/MX
      BINY = (YMAX-YMIN)/MY
      MVN = 0
      MVX = 0.
      MVY = 0.
      MV2 = 0.
      
      OPEN (UNIT=LP_DAT,  STATUS='UNKNOWN', FORM='FORMATTED', &
            FILE="par.dat")
      READ (LP_DAT,*) NP
      DO I=1,NP
      READ (LP_DAT,*) XI,YI,VXI,VYI
      NX = FLOOR((XI-XMIN)/BINX)+1
      NY = FLOOR((YI-YMIN)/BINY)+1
      IF (NX.GE.1.AND.NX.LE.MX.AND.NY.GE.1.AND.NY.LE.MY) THEN 
      MVN(NX,NY) = MVN(NX,NY)+1
      MVX(NX,NY) = MVX(NX,NY)+VXI
      MVY(NX,NY) = MVY(NX,NY)+VYI
      MV2(NX,NY) = MV2(NX,NY)+VXI*VXI+VYI*VYI
      ENDIF
      ENDDO
      CLOSE (LP_DAT)

      OPEN (UNIT=LP_OUT,  STATUS='UNKNOWN', FORM='FORMATTED', &
            FILE="mesh.out")
      OPEN (UNIT=LP_TEMP,  STATUS='UNKNOWN', FORM='FORMATTED', &
            FILE="temp.out")
      OPEN (UNIT=LP_DEN,  STATUS='UNKNOWN', FORM='FORMATTED', &
            FILE="den.out")
      NSUM = 0
      TSUM = 0.
      VXSUM = 0.
      VYSUM = 0.
      DO NY=1,MY
      DO NX=1,MX
      NSUM = NSUM+MVN(NX,NY)
      VXSUM = VXSUM+MVX(NX,NY)
      VYSUM = VYSUM+MVY(NX,NY)
      IF (MVN(NX,NY).GT.0) THEN 
      MVX(NX,NY) = MVX(NX,NY)/MVN(NX,NY)
      MVY(NX,NY) = MVY(NX,NY)/MVN(NX,NY)
      MV2(NX,NY) = MV2(NX,NY)/MVN(NX,NY)
      ENDIF
      MV2(NX,NY) = (MV2(NX,NY)-MVX(NX,NY)**2-MVY(NX,NY)**2)/2.
!     print *, tsum, MV2(NX,NY), MVN(NX,NY)
      TSUM = TSUM+MV2(NX,NY)*MVN(NX,NY)
      WRITE (LP_OUT,'(5F10.4)') XMIN+(NX-0.5)*BINX,YMIN+(NY-0.5)*BINY, &
                                MVX(NX,NY),MVY(NX,NY),MV2(NX,NY)
      ENDDO
      WRITE (LP_DEN,'(100F10.3)') MVN(:,NY)/BINX/BINY
      WRITE (LP_TEMP,'(100F8.5)') MV2(:,NY)
      ENDDO
      CLOSE (LP_OUT)
      CLOSE (LP_DEN)
      CLOSE (LP_TEMP)
      
      PRINT *, 'In a box centered at',(XMIN+XMAX)/2.,(YMIN+YMAX)/2.
      PRINT *, 'with sizex =', XMAX-XMIN, '   sizey =',YMAX-YMIN
      PRINT *, 'There are', NSUM, ' atoms inside.'
      PRINT *, 'Average VX = ',VXSUM/NSUM
      PRINT *, 'Average VY = ',VYSUM/NSUM
      PRINT *, 'Average T = ',TSUM/NSUM
      
      STOP
      END PROGRAM mesh