!     ----------------------------------------------------------------
!     Compile as:  
!     f90 -O3 -n32 -r10000 -OPT:fast_io=ON prop.f90
!     ----------------------------------------------------------------
      PROGRAM property
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: X,Y,VX,VY
      DATA LP_DAT,LP_OUT /25,26/      
      
      XMIN = 0.D0
      XMAX = 150.D0
      YMIN = 0.D0
      YMAX = 50.D0

      NIN  = 0
!     FIRST COUNT HOW MANY ATOMS ARE INSIDE THE REGION:
      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
!     SELECT A REGION
      IF (XI.GT.XMIN.AND.XI.LT.XMAX) THEN   
      NIN = NIN+1
      ENDIF
      ENDDO
      CLOSE (LP_DAT)

      ALLOCATE (X(NIN),Y(NIN),VX(NIN),VY(NIN))

      IP = 0
      VXSUM = 0.
      VYSUM = 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
      IF (XI.GT.XMIN.AND.XI.LT.XMAX) THEN 
      IP = IP+1
      X(IP) = XI 
      Y(IP) = YI 
      VX(IP) = VXI 
      VY(IP) = VYI
      VXSUM = VXSUM+VX(IP)
      VYSUM = VYSUM+VY(IP)
      ENDIF
      ENDDO
      CLOSE (LP_DAT)

      PRINT *, 'In a box centered at',(XMIN+XMAX)/2.,(YMIN+YMAX)/2.
      PRINT *, 'with sizex =', XMAX-XMIN, '   sizey =',YMAX-YMIN
      PRINT *, 'There are', NIN, ' atoms inside.'
      PRINT *, 'Average VX = ',VXSUM/NIN
      PRINT *, 'Average VY = ',VYSUM/NIN

      TMPSUM = 0.
      DO IP=1,NIN
      VX(IP) = VX(IP)-VXSUM/NIN
      VY(IP) = VY(IP)-VYSUM/NIN
      TMPSUM = TMPSUM+VX(IP)*VX(IP)+VY(IP)*VY(IP)
      ENDDO
      PRINT *, 'Average T = ',TMPSUM/NIN/2.
      
      DEALLOCATE (X,Y,VX,VY)
      STOP
      END PROGRAM property