Changeset 1661 in lmdz_wrf

Timestamp:

Oct 4, 2017, 12:27:56 AM (8 years ago)

Author:

lfita

Message:

Finally working version!

• Three different methods have been added to determine the correct coordinate of the trajectory on multiple coincidencies 'mean': get coordinates from the areal-weighted mean of the centers of the given polygons and their areas 'largest': get the coorindates of the largest polygon 'closest': get the coordinates of the closest polygon

Location:

trunk/tools

Files:

• module_scientific.f90 (modified) (22 diffs)
• trajectories_overlap.f90 (modified) (4 diffs)

trunk/tools/module_scientific.f90

@@ -58 +58 @@
 
     INTEGER, INTENT(in)                                  :: funit, dt, itrack
-    REAL(r_k), DIMENSION(4,dt), INTENT(out)              :: ftrack
+    REAL(r_k), DIMENSION(5,dt), INTENT(out)              :: ftrack
 
 ! Local
@@ -79 +79 @@
     DO WHILE (.NOT.found)
 
-      READ(funit,10)ftrack(1,1), Str1, ((ftrack(i,j),Str1,i=2,4),Str1,j=1,dt)
+      READ(funit,10)ftrack(1,1), Str1, ((ftrack(i,j),Str1,i=2,5),Str1,j=1,dt)
       IF (INT(ftrack(1,1)) == itrack) THEN
         ftrack(1,2:dt) = ftrack(1,1)
@@ -92 +92 @@
     RETURN
 
- 10 FORMAT(I10000000,1x,A1,1x,10000000(3(F20.10,A1),A1))
+ 10 FORMAT(I10000000,1x,A1,1x,10000000(4(F20.10,A1),A1))
 
   END SUBROUTINE read_finaltrack_ascii
@@ -102 +102 @@
 
     INTEGER, INTENT(in)                                  :: funit, dt
-    REAL(r_k), DIMENSION(4,dt), INTENT(in)               :: ftrack
+    REAL(r_k), DIMENSION(5,dt), INTENT(in)               :: ftrack
 
 ! Local
@@ -113 +113 @@
 
     fname = 'write_finaltrack_ascii'
-    WRITE(funit,10)INT(ftrack(1,1)), ';', ((ftrack(i,j), ',', i=2,4), ':', j=1,dt)
+    WRITE(funit,10)INT(ftrack(1,1)), ';', ((ftrack(i,j), ',', i=2,5), ':', j=1,dt)
 
     RETURN
 
- 10 FORMAT(I10,1x,A1,1x,10000000(3(F20.10,A1),A1))
+ 10 FORMAT(I10,1x,A1,1x,10000000(4(F20.10,A1),A1))
 
   END SUBROUTINE write_finaltrack_ascii
@@ -373 +373 @@
 
   SUBROUTINE poly_overlap_tracks_area_ascii(dbg, compute, dx, dy, dt, minarea, inNallpolys, allpolys, &
-    ctrpolys, areapolys, Nmaxpoly, Nmaxtracks)
+    ctrpolys, areapolys, Nmaxpoly, Nmaxtracks, methodmulti)
 ! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum 
 !   overlaping/coincidence filtrered by a minimal area writting theoutput on an ASCII file (memory limitations)
@@ -380 +380 @@
 
     LOGICAL, INTENT(in)                                  :: dbg
-    CHARACTER(LEN=*), INTENT(in)                         :: compute
+    CHARACTER(LEN=*), INTENT(in)                         :: compute, methodmulti
     INTEGER, INTENT(in)                                  :: dx, dy, dt, Nmaxpoly, Nmaxtracks
     INTEGER, DIMENSION(dt), INTENT(in)                   :: inNallpolys
@@ -398 +398 @@
     INTEGER, DIMENSION(Nmaxpoly)                         :: prevID, currID
     REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,2)        :: tracks
-    REAL(r_k), DIMENSION(4,dt)                           :: finaltracks
+    REAL(r_k), DIMENSION(5,dt)                           :: finaltracks
     INTEGER, DIMENSION(:), ALLOCATABLE                   :: coins
     INTEGER, DIMENSION(:,:), ALLOCATABLE                 :: coinsNpts
@@ -415 +415 @@
     INTEGER                                              :: idtrack, maxtrack
     REAL(r_k), DIMENSION(5,Nmaxpoly,dt)                  :: singletrack
+    REAL(r_k)                                            :: totArea, dist, mindist, maxarea, areai
 
 !!!!!!! Variables
@@ -428 +429 @@
 ! Nmaxpoly: Maximum possible number of polygons
 ! Nmaxtracks: maximum number of tracks
+! methodmulti: methodology to follow when multiple polygons are given for the same track
+!   'mean': get coordinates from the areal-weighted mean of the centers of the given polygons and their areas
+!   'largest': get the coorindates of the largest polygon
+!   'closest': get the coordinates of the closest polygon
 
     fname = 'poly_overlap_tracks_area_ascii'
@@ -475 +480 @@
     END SELECT
 
+    ! Checking multi-polygon methodology
+    IF ( (TRIM(methodmulti) /= 'mean') .AND. (TRIM(methodmulti) /= 'largest') .AND.                   &
+      (TRIM(methodmulti) /= 'closest')) THEN
+      msg= "methodology for multiple-polygons: '"//TRIM(methodmulti)//"' not ready" // NEW_LINE('a')//&
+        " available ones: 'mean', 'largest', 'closest'"
+      CALL ErrMsg(msg, fname, -1)
+    END IF
+
     IF (dooverlap) THEN
       ! ASCII file for all the polygons and their previous associated one
@@ -585 +598 @@
         msg="Problems allocating 'coinsNpts'"
         CALL ErrMsg(msg,fname,ierr)
-
-        PRINT *,'  Lluis max(searchpolys):', MAXVAL(searchpolys), ' Nsearch:', Nsearch
 
         CALL coincidence_all_polys_area(dbg, dx, dy, Nmeet, currID, meetpolys, Ccurrpolys(:,1:Nmeet),   &
@@ -719 +730 @@
       ! First time-step. Take all polygons
       itrack = 0
-      tracks = 0.
+      tracks = zeroRK
       Ntracks = 0
       it = 1
@@ -776 +787 @@
                 tracks(1,iip,ictrack,it) = itrack
                 tracks(2,iip,ictrack,it) = iiprev
-                ixp = ctrpolys(1,iiprev,iit)
-                iyp = ctrpolys(2,iiprev,iit)
-                tracks(3,iip,ictrack,it) = ixp
-                tracks(4,iip,ictrack,it) = iyp
+                tracks(3,iip,ictrack,it) = ctrpolys(1,iiprev,iit)
+                tracks(4,iip,ictrack,it) = ctrpolys(2,iiprev,iit)
                 tracks(5,iip,ictrack,it) = iit
               END DO
@@ -800 +809 @@
               tracks(1,j,ictrack,it) = maxtrack
               tracks(2,j,ictrack,it) = iiprev
-              ixp = ctrpolys(1,iiprev,iit)
-              iyp = ctrpolys(2,iiprev,iit)
-              tracks(3,j,ictrack,it) = ixp
-              tracks(4,j,ictrack,it) = iyp
+              tracks(3,j,ictrack,it) = ctrpolys(1,iiprev,iit)
+              tracks(4,j,ictrack,it) = ctrpolys(2,iiprev,iit)
               tracks(5,j,ictrack,it) = iit
             END DO
@@ -820 +827 @@
         CALL write_overlap_tracks_ascii(ftrackunit,iit,Nmaxpoly,Ntracks(it),tracks(:,:,1:Ntracks(it),it))
         ! Re-initializing for the next time-step
-        tracks(:,:,:,it-1) = 0
+        tracks(:,:,:,it-1) = zeroRK
         Ntracks(it-1) = Ntracks(it)
         tracks(:,:,1:Ntracks(it-1),it-1) = tracks(:,:,1:Ntracks(it),it)
         Ntracks(it) = 0
-        tracks(:,:,:,it) = 0
+        tracks(:,:,:,it) = zeroRK
 
       END DO timesteps
@@ -849 +856 @@
       CALL ErrMsg(msg,fname,ios)
 
-      finaltracks = 0.
+      finaltracks = zeroRK
 
       DO itt=1, Nmaxtracks
@@ -869 +876 @@
 
         finaltracks = zeroRK
-        finaltracks(1,:) = itrack*1.
+        finaltracks(1,:) = itrack*oneRK
         DO it =1, dt      
           Nprev = COUNT(INT(singletrack(2,:,it)) /= zeroRK)
           IF (Nprev /= 0) THEN
-            finaltracks(2,it) = SUM(singletrack(3,:,it))/Nprev*1.
-            finaltracks(3,it) = SUM(singletrack(4,:,it))/Nprev*1.
-            finaltracks(4,it) = it*1.
+            finaltracks(5,it) = it*oneRK
+            IF (TRIM(methodmulti) == 'largest') THEN
+              maxarea = -10.*oneRK
+              DO ip=1, Nprev
+                IF (areapolys(singletrack(2,ip,it),it) > maxarea) THEN
+                  maxarea = areapolys(singletrack(2,ip,it),it)
+                  i = ip
+                END IF
+              END DO
+              IF (dbg) THEN
+                PRINT *,'  Determine the trajectory coordinates to the largest polygon:', i,          &
+                  ' area:', maxarea
+              END IF
+              finaltracks(2,it) = singletrack(2,i,it)*oneRK
+              finaltracks(3,it) = singletrack(3,i,it)
+              finaltracks(4,it) = singletrack(4,i,it)
+            ELSE IF (TRIM(methodmulti) == 'closest') THEN
+              IF (it > 1) THEN
+                mindist = 10000000.*oneRK
+                DO ip=1, Nprev
+                  dist = SQRT((singletrack(3,ip,it)-finaltracks(3,it-1))**2 +                         &
+                    (singletrack(4,ip,it)-finaltracks(4,it-1))**2 )
+                  IF (dist < mindist) THEN
+                    mindist = dist
+                    i = ip
+                  END IF
+                END DO
+                finaltracks(2,it) = singletrack(3,i,it)*oneRK
+                finaltracks(3,it) = singletrack(3,i,it)
+                finaltracks(4,it) = singletrack(4,i,it)
+                IF (dbg) THEN
+                  PRINT *,'  Determine the trajectory coordinates to the closest previous polygon:',i,&
+                    ' distance:', mindist
+                END IF
+              ELSE
+                maxarea = -10.*oneRK
+                DO ip=1, Nprev
+                  IF (areapolys(singletrack(2,ip,it),it) > maxarea) THEN
+                    maxarea = areapolys(singletrack(2,ip,it),it)
+                    i = ip
+                  END IF
+                END DO
+                IF (dbg) THEN
+                  PRINT *, '  Determine the trajectory coordinates to the largest polygon:', i,        &
+                    ' area:', maxarea, ' at the first time-step then to the closest'
+                END IF
+                finaltracks(2,it) = i*oneRK
+                finaltracks(3,it) = singletrack(3,i,it)
+                finaltracks(4,it) = singletrack(4,i,it)              
+              END IF
+            ELSE
+              totArea = zeroRK
+              finaltracks(2,it) = -oneRK
+              finaltracks(3,it) = zeroRK
+              finaltracks(4,it) = zeroRK
+              DO ip=1, Nprev
+                areai = areapolys(singletrack(2,ip,it),it)
+                totArea = totArea + areai
+                finaltracks(3,it) = finaltracks(3,it) + singletrack(3,ip,it)*areai
+                finaltracks(4,it) = finaltracks(4,it) + singletrack(4,ip,it)*areai
+              END DO
+              finaltracks(3,it) = finaltracks(3,it)/totArea
+              finaltracks(4,it) = finaltracks(4,it)/totArea
+              IF (dbg) THEN
+                PRINT *,'  Determine the trajectory coordinates to the area-averaged polygon ' //     &
+                  ' total area:', totArea
+              END IF
+
+            END IF
+
           END IF
         END DO
@@ -933 +1007 @@
     LOGICAL                                              :: Indeppolychained
     INTEGER                                              :: itrack, ictrack
-    INTEGER                                              :: ixp, iyp, ttrack
+    REAL(r_k)                                            :: ixp, iyp
+    INTEGER                                              :: ttrack
     INTEGER, DIMENSION(dt)                               :: Ntracks
     INTEGER                                              :: idtrack, maxtrack
@@ -1252 +1327 @@
     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
+!   In case of multiple coincidencies, the closest and then the largest is taken filtrered by a minimal area
 !   Here the difference is that the index does not coincide with its ID
 
@@ -1399 +1474 @@
         IF (poly(i,j) .AND. polys(i,j) .NE. 0) THEN
           IF (.NOT.ANY(IDpoly == polys(i,j))) THEN
-            IF ((ip) > 0) PRINT *,'  Lluis ip:', ip, ' Npoly:', Npoly, ' ID:', polys(i,j), ' IDpoly:', IDpoly(1:ip), &
-              ' coinNpts:', coinNpts(ip)
             ip = ip + 1
             IDpoly(ip) = polys(i,j)
@@ -1725 +1798 @@
     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
+!   In case of multiple coincidencies, the closest and then the largest is taken
 
     IMPLICIT NONE

trunk/tools/trajectories_overlap.f90

@@ -21 +21 @@
   CHARACTER(len=50)                                      :: polyfilename
   CHARACTER(len=50)                                      :: poly2Dn, polyNn, polywctrn, polywarean,   &
-    varnresx, varnresy, timen, projectn, waycomp
+    varnresx, varnresy, timen, projectn, waycomp, waymulti
   LOGICAL                                                :: debug
 
@@ -42 +42 @@
 
   NAMELIST /files/ polyfilename, poly2Dn, polyNn, polywctrn, polywarean, varnresx, varnresy, timen
-  NAMELIST /config/ missing_value, minarea, projectn, waycomp, debug
+  NAMELIST /config/ missing_value, minarea, projectn, waycomp, waymulti, debug
 
 !!!!!!!
@@ -63 +63 @@
 !   'scratch': everything from the beginning
 !   'continue': skipt that parts which already have the ascii file written
+! waymulti: methodology to follow when multiple polygons are given for the same track
+!   'mean': get coordinates from the areal-weighted mean of the centers of the given polygons and their areas
+!   'largest': get the coorindates of the largest polygon
+!   'closest': get the coordinates of the closest polygon
 
   fname = 'trajectories_overlap'
@@ -213 +217 @@
   ! Computing trajectories
   CALL poly_overlap_tracks_area_ascii(debug, waycomp, dimx, dimy, dimt, minarea, polyN, polys,        &
-    polywctrs, polyas, Nmaxpoly, Nmaxtrack)
+    polywctrs, polyas, Nmaxpoly, Nmaxtrack, waymulti)
 
   ! Writting trajectory files