Changeset 1652 in lmdz_wrf

Timestamp:

Sep 21, 2017, 7:32:58 PM (8 years ago)

Author:

lfita

Message:

Adding `area'-minimal based system of tracking

File:

• trunk/tools/module_scientific.f90 (modified) (6 diffs)

trunk/tools/module_scientific.f90

@@ -7 +7 @@
 ! coincidence_all_polys: Subtourine to determine which is the coincident polygon when a boolean polygon is provided 
 !   to a map of integer polygons
+! coincidence_all_polys_area: Subtourine to determine which is the coincident polygon when a boolean polygon is provided 
+!   to a map of integer polygons filtrered by a minimal area
 ! coincidence_poly: Subtourine to determine which is the coincident polygon when a boolean polygon is provided 
 !   to a map of integer polygons
+! coincidence_poly_area: Subtourine to determine which is the coincident polygon when a boolean polygon is provided 
+!   to a map of integer polygons filtrered by a minimal area
 ! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
 ! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
@@ -22 +26 @@
 ! poly_overlap_tracks: Subroutine to determine tracks of a series of consecutive 2D field with polygons using 
 !   maximum overlaping/coincidence! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k)
+! poly_overlap_tracks_area: Subroutine to determine tracks of a series of consecutive 2D field with polygons using 
+!   maximum overlaping/coincidence filtrered by a minimal area
 ! quantilesR_K: Subroutine to provide the quantiles of a given set of values of type real 'r_k'
 ! rand_sample: Subroutine to randomly sample a range of indices
@@ -38 +44 @@
 
 
-  SUBROUTINE poly_overlap_tracks(dbg, dx, dy, dt, Nallpolys, allpolys, ctrpolys, areapolys, Nmaxpoly, &
-    Nmaxtracks, tracks, finaltracks)
+  SUBROUTINE poly_overlap_tracks_area(dbg, dx, dy, dt, minarea, inNallpolys, allpolys, ctrpolys,      &
+    areapolys, Nmaxpoly, Nmaxtracks, tracks, finaltracks)
+! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum 
+!   overlaping/coincidence filtrered by a minimal area
+
+    IMPLICIT NONE
+
+    LOGICAL, INTENT(in)                                  :: dbg
+    INTEGER, INTENT(in)                                  :: dx, dy, dt, Nmaxpoly, Nmaxtracks
+    INTEGER, DIMENSION(dt), INTENT(in)                   :: inNallpolys
+    INTEGER, DIMENSION(dx,dy,dt), INTENT(in)             :: allpolys
+    REAL(r_k), INTENT(in)                                :: minarea
+    REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in)      :: ctrpolys
+    REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in)        :: areapolys
+    REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,dt),                                                   &
+      INTENT(out)                                        :: tracks
+    REAL(r_k), DIMENSION(4,Nmaxtracks,dt), INTENT(out)   :: finaltracks
+
+! Local
+    INTEGER                                              :: i, j, ip, it, iip, itt
+    INTEGER                                              :: ierr
+    REAL(r_k), DIMENSION(2,Nmaxpoly,dt)                  :: Actrpolys
+    REAL(r_k), DIMENSION(Nmaxpoly,dt)                    :: Aareapolys
+    INTEGER, DIMENSION(dt)                               :: Nallpolys
+    INTEGER, DIMENSION(dx,dy)                            :: meetpolys, searchpolys
+    INTEGER, DIMENSION(Nmaxpoly)                         :: coincidencies
+    INTEGER, DIMENSION(Nmaxpoly)                         :: prevID, currID
+    INTEGER, DIMENSION(:), ALLOCATABLE                   :: coins
+    INTEGER, DIMENSION(:,:), ALLOCATABLE                 :: coinsNpts
+    INTEGER, DIMENSION(Nmaxpoly,dt)                      :: polycoincidencies
+    INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt)             :: coincidenciesNpts
+    INTEGER                                              :: Nmeet, Nsearch, Nindep
+    INTEGER, DIMENSION(dt)                               :: Nindeppolys
+    CHARACTER(len=5)                                     :: NcoinS
+    INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt)             :: polysIndep
+    INTEGER, DIMENSION(Nmaxpoly,dt)                      :: NpolysIndep
+    INTEGER, DIMENSION(Nmaxpoly,dt)                      :: SpolysIndep
+    INTEGER                                              :: iindep, iiprev
+    INTEGER                                              :: Nprev, NNprev, Ntprev
+    LOGICAL                                              :: Indeppolychained
+    INTEGER                                              :: itrack, ictrack
+    INTEGER                                              :: ixp, iyp, ttrack
+    INTEGER, DIMENSION(dt)                               :: Ntracks
+    INTEGER                                              :: idtrack, maxtrack
+
+!!!!!!! Variables
+! dx,dy,dt: space/time dimensions
+! Nallpolys: Vector with the number of polygons at each time-step
+! allpolys: Series of 2D field with the polygons
+! minarea: minimal area (in same units as areapolys) to perform the tracking 
+! ctrpolys: center of the polygons
+! areapolys: area of the polygons
+! Nmaxpoly: Maximum possible number of polygons
+! Nmaxtracks: maximum number of tracks
+! tracks: series of consecutive polygons
+! trackperiod: period of the track in time-steps
+
+    fname = 'poly_overlap_tracks_area'
+
+    IF (dbg) PRINT *,TRIM(fname)
+
+    polycoincidencies = fillvalI
+    coincidenciesNpts = fillvalI
+    ! Number of independent polygons by time step
+    Nindeppolys = 0
+    ! Number of polygons attached to each independent polygons by time step
+    NpolysIndep = 0
+    ! ID of searching polygon attached to each independent polygons by time step
+    SpolysIndep = 0
+    ! ID of polygons attached to each independent polygons by time step
+    polysIndep = 0
+    ! ID of polygons from previous time-step
+    prevID = 0
+    ! ID of polygons from current time-step
+    currID = 0
+
+    ! First time-step all are independent polygons
+    it = 1
+    Nmeet = inNallpolys(it)
+    Nindeppolys(it) = Nmeet
+    ip = 0
+    meetpolys = allpolys(:,:,it)
+    DO i=1, Nmeet
+      IF (areapolys(i,it) >= minarea) THEN
+        ip = ip + 1
+        SpolysIndep(ip,it) = i
+        currID(ip) = i
+        NpolysIndep(ip,it) = 1
+        polysIndep(ip,1,it) = i
+        Actrpolys(1,ip,it) = ctrpolys(1,i,it)
+        Actrpolys(2,ip,it) = ctrpolys(2,i,it)
+        Aareapolys(ip,it) = areapolys(i,it)
+      ELSE
+        WHERE (meetpolys == i)
+          meetpolys = 0
+        END WHERE
+      END IF
+    END DO
+    Nallpolys(1) = ip
+    Nindeppolys(1) = ip
+
+    ! Starting step
+    it = 0
+    IF (dbg) THEN
+      PRINT *,'  time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',0
+      PRINT *,'    number of independent polygons:', Nindeppolys(it+1)
+      PRINT *,'    indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______'
+      DO i=1, Nindeppolys(it+1)
+        PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':',                                     &
+          polysIndep(i,1:NpolysIndep(i,it+1),it+1)
+      END DO
+    END IF
+
+    ! Looking for the coincidencies at each time step
+    DO it=1, dt-1
+      ! Number of times that a polygon has a coincidence
+      coincidencies = 0
+
+      Nmeet = inNallpolys(it+1)
+      searchpolys = meetpolys
+      meetpolys = allpolys(:,:,it+1)
+      prevID = 0
+      prevID = currID
+      ip = 0
+
+      DO i=1, Nmeet
+        IF (areapolys(i,it+1) >= minarea) THEN
+          ip = ip + 1
+          currID(ip) = i
+          Actrpolys(1,ip,it+1) = ctrpolys(1,i,it+1)
+          Actrpolys(2,ip,it+1) = ctrpolys(2,i,it+1)
+          Aareapolys(ip,it+1) = areapolys(i,it+1)
+        ELSE
+          WHERE (meetpolys == i)
+            meetpolys = 0
+          END WHERE
+        END IF
+      END DO
+      Nallpolys(it+1) = ip
+      Nindeppolys(it+1) = ip
+
+      Nmeet = Nallpolys(it+1)
+      ! Looking throughout the independent polygons
+      Nsearch = Nindeppolys(it)
+
+      IF (ALLOCATED(coins)) DEALLOCATE(coins)
+      ALLOCATE(coins(Nmeet), STAT=ierr)
+      msg="Problems allocating 'coins'"
+      CALL ErrMsg(msg,fname,ierr)
+
+      IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts)
+      ALLOCATE(coinsNpts(Nmeet, Nsearch), STAT=ierr)
+      msg="Problems allocating 'coinsNpts'"
+      CALL ErrMsg(msg,fname,ierr)
+
+      CALL coincidence_all_polys_area(dbg, dx,dy, Nmeet, currID, meetpolys, Actrpolys(:,1:Nmeet,it+1),&
+        Nsearch, prevID, searchpolys, Actrpolys(:,1:Nsearch,it), Aareapolys(1:Nsearch,it), coins,     &
+        coinsNpts)
+
+      polycoincidencies(1:Nmeet,it+1) = coins
+      coincidenciesNpts(1:Nmeet,1:Nsearch,it+1) = coinsNpts
+
+      ! Counting the number of times that a polygon has a coincidency
+      IF (dbg) THEN
+        PRINT *,'  Coincidencies for the given time-step:', it+1,' _______'
+        DO i=1, Nmeet
+          PRINT *,currID(i), coins(i),' N search pts:', coinsNpts(i,:)
+        END DO
+      END IF
+
+      ! Looking for the same equivalencies
+      Nindep = 0
+      DO i=1, Nmeet
+        IF (coins(i) == -1) THEN
+          Nindep = Nindep + 1
+          SpolysIndep(Nindep,it+1) = -1
+          NpolysIndep(Nindep,it+1) = NpolysIndep(Nindep,it+1) + 1
+          polysIndep(Nindep,NpolysIndep(Nindep,it+1),it+1) = currID(i)
+        ELSE IF (coins(i) == -9) THEN
+          WRITE(NcoinS,'(I5)')coins(i)
+          msg="coins= "//TRIM(NcoinS)//" This is an error. One should have always only one " //      &
+            "coincidence of polygon"
+          CALL ErrMsg(msg, fname, -1)
+        ELSE
+          ! Looking for coincidences with previous independent polygons
+          DO ip=1, Nsearch
+            ! Looking into the polygons associated
+            NNprev = NpolysIndep(ip,it)
+            DO j=1, NNprev
+              IF (coins(i) == polysIndep(ip,j,it)) THEN
+                IF (coincidencies(ip) == 0) THEN
+                  Nindep = Nindep + 1
+                  SpolysIndep(Nindep,it+1) = coins(i)
+                END IF
+                ! Which index corresponds to this coincidence?
+                iindep = Index1DArrayI(SpolysIndep(1:Nmeet,it+1), Nindep, coins(i))
+                IF (iindep < 0) THEN
+                  msg = 'Wrong index! There must be an index here'
+                  CALL ErrMsg(msg,fname,iindep)
+                END IF
+                coincidencies(ip) = coincidencies(ip) + 1
+                NpolysIndep(iindep,it+1) = NpolysIndep(iindep,it+1) + 1
+                polysIndep(iindep,NpolysIndep(iindep,it+1),it+1) = currID(i)
+                EXIT
+              END IF
+            END DO
+          END DO
+        END IF
+      END DO
+      Nindeppolys(it+1) = Nindep
+
+      IF (dbg) THEN
+        PRINT *,'  time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',Nsearch
+        PRINT *,'    number of independent polygons:', Nindeppolys(it+1)
+        PRINT *,'    indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______'
+        DO i=1, Nindeppolys(it+1)
+          PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':',                                   &
+            polysIndep(i,1:NpolysIndep(i,it+1),it+1)
+        END DO
+      END IF
+    END DO
+
+    IF (dbg) THEN
+      PRINT *,  'Coincidencies to connect _______'
+      DO it=1, dt
+        PRINT *,'  it:', it, ' Nindep:', Nindeppolys(it)
+        PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs'
+        DO ip=1, Nindeppolys(it)
+          PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':',                         &
+            polysIndep(ip,1:NpolysIndep(ip,it),it)
+        END DO
+      END DO
+
+    END IF
+
+    ! Trajectories
+    ! It should be done following the number of 'independent' polygons
+    ! One would concatenate that independent polygons which share IDs from one step to another
+
+    ! First time-step. Take all polygons
+    itrack = 0
+    tracks = 0.
+    Ntracks = 0
+    DO ip=1, Nindeppolys(1)
+      itrack = itrack + 1
+      tracks(1,1,itrack,1) = itrack*1.
+      tracks(2,1,itrack,1) = SpolysIndep(ip,1)
+      tracks(3,1,itrack,1) = ctrpolys(1,ip,1)
+      tracks(4,1,itrack,1) = ctrpolys(2,ip,1)
+      tracks(5,1,itrack,1) = 1
+      Ntracks(1) = Ntracks(1) + 1
+    END DO
+
+    ! Looping allover already assigned tracks
+    timesteps: DO it=2, dt
+      IF (dbg) PRINT *,'track-timestep:', it, 'N indep polys:', Nindeppolys(it)
+      ! Indep polygons current time-step
+      current_poly: DO i=1, Nindeppolys(it)
+        IF (dbg) PRINT *,'  curent poly:', i, 'Prev poly:', SpolysIndep(i,it), ' N ass. polygons:',   &
+          NpolysIndep(i,it), 'ass. poly:', polysIndep(i,1:NpolysIndep(i,it),it)
+        Indeppolychained = .FALSE.
+
+        ! Number of tracks previous time-step
+        ! Looping overall
+        it1_tracks: DO itt=1, Ntracks(it-1)
+          itrack = tracks(1,1,itt,it-1)
+          ! Number polygons ID assigned
+          Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0)
+          IF (dbg) PRINT *,itt,'  track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1)
+
+          ! Looking for coincidencies
+          DO iip=1, Ntprev
+            IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN
+              Indeppolychained = .TRUE.
+              IF (dbg) PRINT *,'    coincidence found by polygon:', tracks(2,iip,itt,it-1)
+              EXIT
+            END IF
+          END DO
+          IF (Indeppolychained) THEN
+            Ntracks(it) = Ntracks(it) + 1
+            ictrack = Ntracks(it)
+            ! Assigning all the IDs to the next step of the track
+            DO iip=1, NpolysIndep(i,it)
+              iiprev = polysIndep(i,iip,it)
+              tracks(1,iip,ictrack,it) = itrack
+              tracks(2,iip,ictrack,it) = iiprev
+              ixp = ctrpolys(1,iiprev,it)
+              iyp = ctrpolys(2,iiprev,it)
+              tracks(3,iip,ictrack,it) = ixp
+              tracks(4,iip,ictrack,it) = iyp
+              tracks(5,iip,ictrack,it) = it
+            END DO
+            EXIT
+          END IF
+          IF (Indeppolychained) EXIT
+        END DO it1_tracks
+
+        ! Creation of a new track
+        IF (.NOT.Indeppolychained) THEN
+          Ntracks(it) = Ntracks(it) + 1
+          ictrack = Ntracks(it)
+          ! ID of new track
+          maxtrack = INT(MAXVAL(tracks(1,:,:,:)*1.))
+          IF (dbg) PRINT *,'  New track!', maxtrack+1
+
+          ! Assigning all the IDs to the next step of the track
+          DO j=1, NpolysIndep(i,it)
+            iiprev = polysIndep(i,j,it)
+            tracks(1,j,ictrack,it) = maxtrack+1
+            tracks(2,j,ictrack,it) = iiprev
+            ixp = ctrpolys(1,iiprev,it)
+            iyp = ctrpolys(2,iiprev,it)
+            tracks(3,j,ictrack,it) = ixp
+            tracks(4,j,ictrack,it) = iyp
+            tracks(5,j,ictrack,it) = it
+          END DO
+        END IF
+
+      END DO current_poly
+
+      IF (dbg) THEN
+        PRINT *,'  At this time-step:', it, ' N trajectories:', Ntracks(it)
+        DO i=1, Ntracks(it)
+          Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0)
+          PRINT *,i ,'ID tracks:', tracks(1,1,i,it), 'ID polygon:', tracks(2,1:Nprev,i,it)
+        END DO
+      END IF
+
+    END DO timesteps
+
+    ! Summarizing trajectories
+    ! When multiple polygons are available, the mean of their central positions determines the position
+
+    finaltracks = 0.
+    maxtrack = MAXVAL(tracks(1,:,:,:))
+
+    DO it=1, dt
+      DO itt=1, Ntracks(it)
+        itrack = INT(tracks(1,1,itt,it))
+        Nprev = COUNT(INT(tracks(2,:,itt,it)) /= 0)
+        PRINT *,'it:', it,'itrack:', itrack, 'Nprev:', Nprev
+        PRINT *,'  coordinates ix iy _______'
+        DO i=1,Nprev
+          PRINT *,'    ',tracks(3,i,itt,it),tracks(4,i,itt,it)
+        END DO
+        finaltracks(1,itrack,it) = itrack*1.
+        finaltracks(2,itrack,it) = SUM(tracks(3,:,itt,it))/Nprev
+        finaltracks(3,itrack,it) = SUM(tracks(4,:,itt,it))/Nprev
+        finaltracks(4,itrack,it) = it*1.
+        PRINT *,'  finaltrack:', finaltracks(:,itrack,it)
+      END DO
+    END DO
+
+    RETURN
+
+  END SUBROUTINE poly_overlap_tracks_area
+
+  SUBROUTINE coincidence_all_polys_area(dbg, dx, dy, Nallpoly, IDallpoly, allpoly, icpolys, Npoly, &
+    IDpolys, polys, cpolys, apolys, polycoins, coinNptss)
+! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
+!   In case of multiple coincidencies, the closest and then the biggest is taken filtrered by a minimal area
+!   Here the difference is that the index does not coincide with its ID
+
+    IMPLICIT NONE
+
+    LOGICAL, INTENT(in)                                  :: dbg
+    INTEGER, INTENT(in)                                  :: dx, dy, Nallpoly, Npoly
+    INTEGER, DIMENSION(dx,dy), INTENT(in)                :: allpoly, polys
+    INTEGER, DIMENSION(Nallpoly), INTENT(in)             :: IDallpoly
+    INTEGER, DIMENSION(Npoly), INTENT(in)                :: IDpolys
+    REAL(r_k), DIMENSION(2,Nallpoly), INTENT(in)         :: icpolys
+    REAL(r_k), DIMENSION(2,Npoly), INTENT(in)            :: cpolys
+    REAL(r_k), DIMENSION(Npoly), INTENT(in)              :: apolys
+    INTEGER, DIMENSION(Nallpoly), INTENT(out)            :: polycoins
+    INTEGER, DIMENSION(Nallpoly,Npoly), INTENT(out)      :: coinNptss
+
+! Local
+    INTEGER                                              :: i, j, ip
+    INTEGER                                              :: maxcorr
+    INTEGER                                              :: Nmaxcorr
+    LOGICAL, DIMENSION(dx,dy)                            :: boolpoly
+    INTEGER                                              :: maxcoin
+    REAL                                                 :: dist, maxcoindist, maxcoinarea
+    INTEGER, DIMENSION(Npoly)                            :: IDcoins
+
+!!!!!!! Variables
+! dx,dy: dimension of the space
+! Nallpoly: Number of polygons to find coincidence
+! allpoly: space with the polygons to meet
+! icpolys: center of the polygons to look for the coincidence
+! Npoly: number of polygons on the 2D space
+! polys: 2D field of polygons identified by their integer number (0 for no polygon)
+! cpolys: center of the polygons
+! apolys: area of the polygons
+! polycoins: coincident polyogn
+!          -1: no-coincidence
+!   1 < Npoly: single coincidence with a given polygon
+!          -9: coincidence with more than one polygon
+! coinNptss: number of points coincident with each polygon
+
+    fname = 'coincidence_all_polys_area'
+    IF (dbg) PRINT *,TRIM(fname)
+
+    DO ip=1, Nallpoly
+      boolpoly = allpoly == IDallpoly(ip)
+      CALL coincidence_poly_area(dbg, dx, dy, boolpoly, Npoly, polys, polycoins(ip), IDcoins,         &
+        coinNptss(ip,:))
+      IF (dbg) PRINT *,'  polygon', IDallpoly(ip), ' coincidence with:', polycoins(ip), 'IDpolys:', IDpolys(1:Npoly)
+
+      ! Coincidence with more than one polygon
+      IF (polycoins(ip) == -9) THEN
+        maxcoindist = -10.
+        maxcoinarea = -10.
+        maxcoin = MAXVAL(coinNptss(ip,:))
+        DO j=1, Npoly
+          IF (coinNptss(ip,j) == maxcoin) THEN
+            dist = SQRT( (icpolys(1,ip)*1.-cpolys(1,j)*1.)**2 + (icpolys(2,ip)*1.-cpolys(2,j)*1.)**2 )
+            IF ( dist > maxcoindist) THEN
+              maxcoindist = dist
+              maxcoinarea = apolys(j)
+              polycoins(ip) = IDcoins(j)
+            ELSE IF ( dist == maxcoindist) THEN
+              IF (apolys(j) > maxcoinarea) THEN
+                polycoins(ip) = IDcoins(j)
+                maxcoinarea = apolys(j)
+              END IF
+            END IF 
+          END IF
+        END DO
+      END IF
+    END DO
+
+    RETURN
+
+  END SUBROUTINE coincidence_all_polys_area
+
+  SUBROUTINE coincidence_poly_area(dbg, dx, dy, poly, Npoly, polys, polycoin, IDpoly, coinNpts)
+! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
+!   Here the difference is that the index does not coincide with its ID
+
+    IMPLICIT NONE
+
+    LOGICAL, INTENT(in)                                  :: dbg
+    INTEGER, INTENT(in)                                  :: dx, dy, Npoly
+    LOGICAL, DIMENSION(dx,dy), INTENT(in)                :: poly
+    INTEGER, DIMENSION(dx,dy), INTENT(in)                :: polys
+    INTEGER, INTENT(out)                                 :: polycoin
+    INTEGER, DIMENSION(Npoly), INTENT(out)               :: IDpoly, coinNpts
+
+! Local
+    INTEGER                                              :: i, j, ip
+    INTEGER                                              :: maxcorr
+    INTEGER                                              :: Nmaxcorr
+
+!!!!!!! Variables
+! dx,dy: dimension of the space
+! poly: bolean polygon to meet
+! Npoly: number of polygons on the 2D space
+! polys: 2D field of polygons identified by their integer number (0 for no polygon)
+! polycoin: coincident polyogn
+!          -1: no-coincidence
+!   1 < Npoly: single coincidence with a given polygon
+!          -9: coincidence with more than one polygon
+! IDpoly: ID of the found polygon
+! coinNpts: number of points coincident with each polygon
+
+    fname = 'coincidence_poly_area'
+    IF (dbg) PRINT *,TRIM(fname)
+
+    IF (dbg) THEN
+      PRINT *,'  Boolean polygon to search coincidences ...'
+      DO i=1,dx
+        PRINT *,poly(i,:)
+      END DO
+
+      PRINT *,'  2D polygons space ...'
+      DO i=1,dx
+        PRINT '(1000(I3,1x))',polys(i,:)
+      END DO
+    END IF
+
+    ! Looking for coincient points for the polygon
+    coinNpts = 0
+    IDpoly = 0
+    ip = 0
+    DO i=1,dx
+      DO j=1,dy
+        IF (poly(i,j) .AND. polys(i,j) .NE. 0) THEN
+          IF (.NOT.ANY(IDpoly == polys(i,j))) THEN
+            ip = ip + 1
+            IDpoly(ip) = polys(i,j)
+          ELSE
+            ip = Index1DarrayI(IDpoly, Npoly, polys(i,j))
+          END IF
+          coinNpts(ip) = coinNpts(ip) + 1
+        END IF
+      END DO
+    END DO
+
+    maxcorr = 0
+    maxcorr = INT(MAXVAL(coinNpts*1.))
+
+    IF (dbg) PRINT *,'  Maximum coincidence:', maxcorr
+    IF (maxcorr == 0) THEN
+      polycoin = -1
+    ELSE
+      Nmaxcorr = 0
+      DO ip=1, Npoly
+        IF (coinNpts(ip) == maxcorr) THEN
+          Nmaxcorr = Nmaxcorr+1
+          polycoin = IDpoly(ip)
+        END IF
+      END DO
+      IF (Nmaxcorr > 1) polycoin = -9
+    END IF
+
+    IF (dbg) THEN
+      PRINT *,'  Coincidences for each polygon _______', Npoly
+      DO ip=1, Npoly
+        PRINT *,'  ',ip, ' ID:', IDpoly(ip),': ', coinNpts(ip)
+      END DO
+    END IF
+
+    RETURN
+
+END SUBROUTINE coincidence_poly_area
+
+  SUBROUTINE poly_overlap_tracks(dbg, dx, dy, dt, minarea, Nallpolys, allpolys, ctrpolys,        &
+    areapolys, Nmaxpoly, Nmaxtracks, tracks, finaltracks)
 ! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum overlaping/coincidence
 
@@ -48 +581 @@
     INTEGER, DIMENSION(dt), INTENT(in)                   :: Nallpolys
     INTEGER, DIMENSION(dx,dy,dt), INTENT(in)             :: allpolys
+    REAL(r_k), INTENT(in)                                :: minarea
     REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in)      :: ctrpolys
     REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in)        :: areapolys
@@ -80 +614 @@
 ! Nallpolys: Vector with the number of polygons at each time-step
 ! allpolys: Series of 2D field with the polygons
+! minarea: minimal area (in same units as areapolys) to perform the tracking 
 ! ctrpolys: center of the polygons
 ! areapolys: area of the polygons
@@ -286 +821 @@
       END DO current_poly
 
-      IF (dbg) PRINT *,'  At this time-step:', it, ' N trajectories:', Ntracks(it)
-      DO i=1, Ntracks(it)
-        Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0)
-        PRINT *,i,'tracks:', tracks(2,1:Nprev,i,it)
-      END DO
+      IF (dbg) THEN
+        PRINT *,'  At this time-step:', it, ' N trajectories:', Ntracks(it)
+        DO i=1, Ntracks(it)
+          Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0)
+          PRINT *,i,'tracks:', tracks(2,1:Nprev,i,it)
+        END DO
+      END IF
 
     END DO timesteps