Index: /trunk/tools/module_scientific.f90
===================================================================
--- /trunk/tools/module_scientific.f90	(revision 1649)
+++ /trunk/tools/module_scientific.f90	(revision 1650)
@@ -5,4 +5,8 @@
 ! all_polygons_properties: Subroutine to determine the properties of all polygons in a 2D field:
 ! borders_matrixL: Subroutine to provide the borders of a logical array (interested in .TRUE.)
+! coincidence_all_polys: Subtourine to determine which is the coincident polygon when a boolean polygon is provided 
+!   to a map of integer polygons
+! coincidence_poly: Subtourine to determine which is the coincident polygon when a boolean polygon is provided 
+!   to a map of integer polygons
 ! 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.
@@ -16,5 +20,6 @@
 ! polygons: Subroutine to search the polygons of a border field. FORTRAN based. 1st = 1!
 ! polygons_t: Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!
-! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k)
+! 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)
 ! 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
@@ -31,4 +36,450 @@
 
   CONTAINS
+
+  SUBROUTINE poly_overlap_tracks(dbg, dx, dy, dt, 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
+
+    IMPLICIT NONE
+
+    LOGICAL, INTENT(in)                                  :: dbg
+    INTEGER, INTENT(in)                                  :: dx, dy, dt, Nmaxpoly, Nmaxtracks
+    INTEGER, DIMENSION(dt), INTENT(in)                   :: Nallpolys
+    INTEGER, DIMENSION(dt,dx,dy), INTENT(in)             :: allpolys
+    REAL(r_k), DIMENSION(dt,Nmaxpoly,2), INTENT(in)      :: ctrpolys
+    REAL(r_k), DIMENSION(dt,Nmaxpoly), INTENT(in)        :: areapolys
+    REAL(r_k), DIMENSION(dt,Nmaxtracks,Nmaxpoly,5),                                                   &
+      INTENT(out)                                        :: tracks
+    REAL(r_k), DIMENSION(dt,Nmaxtracks,4), INTENT(out)   :: finaltracks
+
+! Local
+    INTEGER                                              :: i, j, ip, it, iip, itt
+    INTEGER                                              :: ierr
+    INTEGER, DIMENSION(dt,Nmaxpoly)                      :: coincidencies, NOcoincidencies,           &
+      MULTIcoincidencies
+    INTEGER, DIMENSION(:), ALLOCATABLE                   :: coins
+    INTEGER, DIMENSION(:,:), ALLOCATABLE                 :: coinsNpts
+    INTEGER, DIMENSION(dt,Nmaxpoly)                      :: polycoincidencies
+    INTEGER, DIMENSION(dt,Nmaxpoly,Nmaxpoly)             :: coincidenciesNpts
+    INTEGER                                              :: Nmeet, Nsearch, Nindep
+    INTEGER, DIMENSION(dt)                               :: Nindeppolys
+    CHARACTER(len=5)                                     :: NcoinS
+    INTEGER, DIMENSION(dt,Nmaxpoly,Nmaxpoly)             :: polysIndep
+    INTEGER, DIMENSION(dt,Nmaxpoly)                      :: NpolysIndep
+    INTEGER, DIMENSION(dt,Nmaxpoly)                      :: SpolysIndep
+    INTEGER                                              :: iindep, iiprev
+    INTEGER, DIMENSION(dt,Nmaxpoly,Nmaxpoly)             :: polysIndepchains
+    INTEGER                                              :: Nprev, polyprev, Nnotchained, Nprevchained
+    INTEGER                                              :: NNprev, Ntprev
+    LOGICAL                                              :: Indeppolychained, prevchained
+    INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly)                :: polysIndepNotchained
+    INTEGER, DIMENSION(2)                                :: trackperiod
+    INTEGER                                              :: itrack, ictrack
+    INTEGER                                              :: ixp, iyp, ttrack
+    INTEGER                                              :: Nnew
+    LOGICAL, DIMENSION(Nmaxpoly)                         :: chained
+    INTEGER, DIMENSION(dt)                               :: Ntracks
+    INTEGER                                              :: idtrack, maxtrack
+    INTEGER, DIMENSION(dt,Nmaxpoly)                      :: tracks_it
+
+!!!!!!! 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
+! 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'
+
+    IF (dbg) PRINT *,TRIM(fname)
+
+    polycoincidencies = fillvalI
+    coincidenciesNpts = fillvalI
+    ! Number of times that a polygon has a coincidence
+    coincidencies = 0
+    ! Polygons without a coincidence
+    NOcoincidencies = 0
+    ! Polygons with multiple coincidences
+    MULTIcoincidencies = 0
+    ! 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
+
+    ! First time-step all are independent polygons
+    it = 1
+    Nmeet = Nallpolys(it)
+    Nindeppolys(it) = Nmeet
+    DO i=1, Nmeet
+      SpolysIndep(it,i) = i
+      NpolysIndep(it,1:Nmeet) = 1
+      polysIndep(it,i,1) = i
+    END DO
+
+    ! Looking for the coincidencies at each time step
+    DO it=1, dt-1
+      Nmeet = Nallpolys(it+1)
+      Nsearch = Nallpolys(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(dbg, dx, dy, Nmeet, allpolys(it+1,:,:), ctrpolys(it+1,1:Nmeet,:),    &
+        Nsearch, allpolys(it,:,:), ctrpolys(it,1:Nsearch,:), areapolys(it,1:Nsearch), coins, coinsNpts)
+
+      polycoincidencies(it+1,1:Nmeet) = coins
+      coincidenciesNpts(it+1,1:Nmeet,1:Nsearch) = 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 *,coins(i),' N search pts:', coinsNpts(i,:)
+        END DO
+      END IF
+
+      Nindep = 0
+      DO i=1, Nmeet
+        IF (coins(i) == -1) THEN
+          Nindep = Nindep + 1
+          NOcoincidencies(it+1,i) = 1
+          SpolysIndep(it+1,Nindep) = -1
+          NpolysIndep(it+1,Nindep) = NpolysIndep(it+1,Nindep) + 1
+          polysIndep(it+1,Nindep,NpolysIndep(it+1,Nindep)) = 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
+          DO ip=1, Nsearch
+            IF (coins(i) == ip) THEN
+              IF (coincidencies(it+1,ip) == 0) THEN
+                Nindep = Nindep + 1
+                SpolysIndep(it+1,Nindep) = ip
+              END IF
+              coincidencies(it+1,ip) = coincidencies(it+1,ip) + 1
+              DO iindep=1, Nindep
+                IF (SpolysIndep(it+1,iindep) == coins(i)) THEN
+                  NpolysIndep(it+1,iindep) = NpolysIndep(it+1,iindep) + 1
+                  polysIndep(it+1,iindep,NpolysIndep(it+1,iindep)) = i
+                END IF
+              END DO
+            END IF
+          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(it+1,i), NpolysIndep(it+1,i), ':',                                   &
+            polysIndep(it+1,i,1:NpolysIndep(it+1,i))
+        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(it,ip), ':',                         &
+            polysIndep(it,ip,1:NpolysIndep(it,ip))
+        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,itrack,1,1) = itrack*1.
+      tracks(1,itrack,1,2) = SpolysIndep(1,ip)
+      tracks(1,itrack,1,3) = ctrpolys(1,ip,1)
+      tracks(1,itrack,1,4) = ctrpolys(1,ip,2)
+      tracks(1,itrack,1,5) = 1
+      tracks_it(1,itrack) = itrack
+      Ntracks(1) = Ntracks(1) + 1
+    END DO
+
+    ! Looping allover already assigned tracks
+    timesteps: DO it=2, dt
+      IF (dbg) PRINT *,'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(it,i), ' N ass. polygons:',   &
+          NpolysIndep(it,i), 'ass. poly:', polysIndep(it,i,1:NpolysIndep(it,i))
+        Indeppolychained = .FALSE.
+
+        ! Number of tracks previous time-step
+        ! Looping overall
+        it1_tracks: DO itt=1, Ntracks(it-1)
+          itrack = tracks(it-1,itt,1,1)
+          ! Number polygons ID assigned
+          Ntprev = COUNT(tracks(it-1,itt,:,2) /= 0)
+          IF (dbg) PRINT *,itt,'  track:', itrack, 'assigned:', tracks(it-1,itt,1:Ntprev,2)
+
+          ! Looking for coincidencies
+          DO iip=1, Ntprev
+            IF (tracks(it-1,itt,iip,2) == SpolysIndep(it,i)) THEN
+              Indeppolychained = .TRUE.
+              IF (dbg) PRINT *,'    coincidence found by polygon:', tracks(it-1,itt,iip,2)
+              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(it,i)
+              iiprev = polysIndep(it,i,iip)
+              tracks(it,ictrack,iip,1) = itrack
+              tracks(it,ictrack,iip,2) = iiprev
+              ixp = ctrpolys(it,iiprev,1)
+              iyp = ctrpolys(it,iiprev,2)
+              tracks(it,ictrack,iip,3) = ixp
+              tracks(it,ictrack,iip,4) = iyp
+              tracks(it,ictrack,iip,5) = it
+            END DO
+            tracks_it(it,ictrack) = itrack
+            EXIT
+          END IF
+        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(it,i)
+            iiprev = polysIndep(it,i,j)
+            tracks(it,ictrack,j,1) = maxtrack+1
+            tracks(it,ictrack,j,2) = iiprev
+            ixp = ctrpolys(it,iiprev,1)
+            iyp = ctrpolys(it,iiprev,2)
+            tracks(it,ictrack,j,3) = ixp
+            tracks(it,ictrack,j,4) = iyp
+            tracks(it,ictrack,j,5) = it
+          END DO
+          tracks_it(it,Ntracks(it)) = maxtrack+1
+        END IF
+
+      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(it,i,:,2)) /= 0)
+        PRINT *,i,'tracks:', tracks(it,i,1:Nprev,2)
+      END DO
+
+    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(it,itt,1,1))
+        Nprev = COUNT(INT(tracks(it,itt,:,2)) /= 0)
+        PRINT *,'it:', it,'itrack:', itrack, 'Nprev:', Nprev
+        finaltracks(it,itrack,1) = itrack*1.
+        finaltracks(it,itrack,2) = SUM(tracks(it,itt,:,3))/Nprev
+        finaltracks(it,itrack,3) = SUM(tracks(it,itt,:,4))/Nprev
+        finaltracks(it,itrack,4) = it*1.
+        PRINT *,'  finaltrack:', finaltracks(it,itrack,:)
+      END DO
+    END DO
+
+    RETURN
+
+  END SUBROUTINE poly_overlap_tracks
+
+  SUBROUTINE coincidence_all_polys(dbg, dx, dy, Nallpoly, allpoly, icpolys, Npoly, 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
+
+    IMPLICIT NONE
+
+    LOGICAL, INTENT(in)                                  :: dbg
+    INTEGER, INTENT(in)                                  :: dx, dy, Nallpoly, Npoly
+    INTEGER, DIMENSION(dx,dy), INTENT(in)                :: allpoly, polys
+    REAL(r_k), DIMENSION(Nallpoly,2), INTENT(in)         :: icpolys
+    REAL(r_k), DIMENSION(Npoly,2), 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
+
+!!!!!!! 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'
+    IF (dbg) PRINT *,TRIM(fname)
+
+    DO ip=1, Nallpoly
+      boolpoly = allpoly == ip
+      CALL coincidence_poly(dbg, dx, dy, boolpoly, Npoly, polys, polycoins(ip), coinNptss(ip,:))
+      IF (dbg) PRINT *,'  polygon', ip, ' coincidence with:', polycoins(ip)
+
+      ! 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(ip,1)*1.-cpolys(j,1)*1.)**2 + (icpolys(ip,2)*1.-cpolys(j,2)*1.)**2 )
+            IF ( dist > maxcoindist) THEN
+              maxcoindist = dist
+              maxcoinarea = apolys(j)
+              polycoins(ip) = j
+            ELSE IF ( dist == maxcoindist) THEN
+              IF (apolys(j) > maxcoinarea) THEN
+                polycoins(ip) = j
+                maxcoinarea = apolys(j)
+              END IF
+            END IF 
+          END IF
+        END DO
+      END IF
+    END DO
+
+    RETURN
+
+  END SUBROUTINE coincidence_all_polys
+
+  SUBROUTINE coincidence_poly(dbg, dx, dy, poly, Npoly, polys, polycoin, coinNpts)
+! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
+
+    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)               :: 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
+! coinNpts: number of points coincident with each polygon
+
+    fname = 'coincidence_poly'
+    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
+    DO i=1,dx
+      DO j=1,dy
+        IF (poly(i,j) .AND. polys(i,j) .NE. 0) coinNpts(polys(i,j)) = coinNpts(polys(i,j)) + 1
+      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 = ip
+        END IF
+      END DO
+      IF (Nmaxcorr > 1) polycoin = -9
+    END IF
+
+    IF (dbg) THEN
+      PRINT *,'  Coincidences for each polygon _______'
+      DO ip=1, Npoly
+        PRINT *,'  ', ip,': ', coinNpts(ip)
+      END DO
+    END IF
+
+    RETURN
+
+END SUBROUTINE coincidence_poly
+
 
   SUBROUTINE all_polygons_properties(dbg, dx, dy, Npoly, polys, lon, lat, values, xres, yres, projN,  &
