Changeset 2340 in lmdz_wrf for trunk/tools

Timestamp:

Feb 15, 2019, 5:50:02 PM (6 years ago)

Author:

lfita

Message:

Adding:

• `continguos_homogene_zones': Subroutine to look for contiguous zones by looking by continuous grid points
• `get_xyconlimits': Subroutine for getting the limits of contiguous values from a given point in a 2D matrix
• `Index1DArrayL': Function to provide the first index of a given value inside a 1D boolean array
• `multi_Index1DArrayL': Subroutine to provide the indices of a given value inside a 1D boolean array

File:

• trunk/tools/module_generic.f90 (modified) (3 diffs)

trunk/tools/module_generic.f90

@@ -3 +3 @@
 
 !!!!!!! Subroutines/Functions
+! continguos_homogene_zones: Subroutine to look for contiguous zones by looking by continuous grid points
 ! freeunit: provides the number of a free unit in which open a file
 ! GetInNamelist: Subroutine to get a paramter from a namelistfile
+! get_xyconlimits: Subroutine for getting the limits of contiguous values from a given point in a 2D matrix
 ! index_list_coordsI: Function to provide the index of a given coordinate within a list of integer coordinates
 ! Index1DArrayI: Function to provide the first index of a given value inside a 1D integer array
 ! Index1DArrayR: Function to provide the first index of a given value inside a 1D real array
 ! Index1DArrayR_K: Function to provide the first index of a given value inside a 1D real(r_k) array
+! Index1DArrayL: Function to provide the first index of a given value inside a 1D boolean array
 ! Index2DArrayR: Function to provide the first index of a given value inside a 2D real array
 ! Index2DArrayR_K: Function to provide the first index of a given value inside a 2D real(r_k) array
 ! Nvalues_2DArrayI: Number of different values of a 2D integer array
 ! mat2DPosition: Function to provide the i, j indices of a given value inside a 2D matrix
+! multi_Index1DArrayL: Subroutine to provide the indices of a given value inside a 1D boolean array
 ! numberTimes: Function to provide the number of times that a given set of characters happen within a string
 ! RangeI: Function to provide a range of d1 values from 'iniv' to 'endv', of integer values in a vector
@@ -174 +178 @@
 
   END FUNCTION Index1DArrayR_K
+
+  INTEGER FUNCTION Index1DArrayL(array1D, d1, val)
+! Function to provide the first index of a given value inside a 1D boolean array
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1
+    LOGICAL, INTENT(in)                                  :: val
+    LOGICAL, DIMENSION(d1), INTENT(in)                   :: array1D
+
+! Local
+    INTEGER                                              :: i
+    CHARACTER(LEN=50)                                    :: fname
+
+    fname = 'Index1DArrayL'
+
+    Index1DArrayL = -1
+
+    DO i=1,d1
+      IF (array1d(i) .EQV. val) THEN
+        Index1DArrayL = i
+        EXIT
+      END IF
+    END DO
+
+  END FUNCTION Index1DArrayL
+
+  SUBROUTINE multi_Index1DArrayL(array1D, d1, val, Ntimes, pos)
+! Subroutine to provide the indices of a given value inside a 1D boolean array
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1
+    LOGICAL, INTENT(in)                                  :: val
+    LOGICAL, DIMENSION(d1), INTENT(in)                   :: array1D
+    INTEGER, INTENT(out)                                 :: Ntimes
+    INTEGER, DIMENSION(d1), INTENT(out)                  :: pos
+
+! Local
+    INTEGER                                              :: i
+    CHARACTER(LEN=50)                                    :: fname
+
+!!!!!!! Variables
+! array1D: 1D Array of values
+! d1: length of array
+! val: Value to look for
+! Ntimes: Number of times val is found within array1D
+! pos: positions of the values of val
+
+    fname = 'multi_Index1DArrayL'
+
+    Ntimes = 0
+    pos = -1
+
+    DO i=1,d1
+      IF (array1d(i) .EQV. val) THEN
+        Ntimes = Ntimes + 1
+        pos(Ntimes) = i
+      END IF
+    END DO
+
+  END SUBROUTINE multi_Index1DArrayL
 
   FUNCTION Index2DArrayR(array2D, d1, d2, val)
@@ -704 +770 @@
   END FUNCTION vectorR_S
 
+  SUBROUTINE continguos_homogene_zones(dx, dy, matvals, Nzones, contzones)
+  ! Subroutine to look for contiguous zones by looking by continuous grid points
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: dx, dy
+    INTEGER, DIMENSION(dx,dy), INTENT(in)                :: matvals
+    INTEGER, INTENT(out)                                 :: Nzones
+    INTEGER, DIMENSION(dx,dy), INTENT(out)               :: contzones
+! Local
+    INTEGER                                              :: i,j, k
+    INTEGER                                              :: ii, ei, ij, ej
+    INTEGER                                              :: Ncont, Nassigned, Ncontmin
+    LOGICAL, DIMENSION(dx,dy)                            :: assigned, notdone
+    INTEGER, DIMENSION(:), ALLOCATABLE                   :: pzones, allzones
+    LOGICAL, DIMENSION(:), ALLOCATABLE                   :: passigns
+
+!!!!!!! Variables
+! dx, dy: shape of the matrix
+! matvals: matrix with the values
+! contzones: homogeneous zones found
+
+    fname = 'continguos_homogene_zones'
+
+    ! Vector to keep track of all zone values
+    IF (ALLOCATED(allzones)) DEALLOCATE(allzones)
+    ALLOCATE(allzones(dx*dy/4))
+    allzones = 0
+
+    assigned = .FALSE.
+    notdone = .TRUE.
+    contzones = -1
+    Nzones = 0
+    DO i=1, dx
+      DO j=1, dy
+        ! First
+        IF (matvals(i,j) /= 0 .AND. notdone(i,j)) THEN
+          CALL get_xyconlimits(dx, dy, matvals, 0, i, j, ii, ei, ij, ej)
+
+          ! Has any point of the rays already been assigned?
+          ! Along x
+          IF (ALLOCATED(pzones)) DEALLOCATE(pzones)
+          ALLOCATE(pzones(ei-ii+1))
+          IF (ALLOCATED(passigns)) DEALLOCATE(passigns)
+          ALLOCATE(passigns(ei-ii+1))
+          passigns = assigned(ii:ei,j)
+          CALL multi_Index1DArrayL(passigns, ei-ii+1, .TRUE., Nassigned, pzones)
+          IF (Nassigned /= 0) THEN
+            Ncontmin = 10000000
+            DO k=1, Nassigned
+              IF (contzones(ii+pzones(k)-1,j) < Ncontmin) Ncontmin = contzones(ii+pzones(k)-1,j)
+            END DO
+            ! If there is more than one assigned value change all values to the minimum one
+            DO k=1, Nassigned
+              IF (contzones(ii+pzones(k)-1,j) /= Ncontmin) THEN
+                WHERE (contzones == contzones(ii+pzones(k)-1,j) .AND. contzones /= Ncontmin)
+                  contzones = Ncontmin
+                END WHERE
+                allzones(contzones(ii+pzones(k)-1,j)) = 0
+              END IF
+            END DO
+            Ncont = Ncontmin
+          END IF
+          ! Along y
+          IF (ALLOCATED(pzones)) DEALLOCATE(pzones)
+          ALLOCATE(pzones(ej-ij+1))
+          IF (ALLOCATED(passigns)) DEALLOCATE(passigns)
+          ALLOCATE(passigns(ej-ij+1))
+          passigns = assigned(i,ij:ej)
+          CALL multi_Index1DArrayL(passigns, ej-ij+1, .TRUE., Nassigned, pzones)
+          IF (Nassigned /= 0) THEN
+            Ncontmin = 10000000
+            DO k=1, Nassigned
+              IF (contzones(i,ij+pzones(k)-1) < Ncontmin) Ncontmin = contzones(i,ij+pzones(k)-1)
+            END DO
+            ! If there is more than one assigned value change all values to the minimum one
+            DO k=1, Nassigned
+              IF (contzones(i,ij+pzones(k)-1) /= Ncontmin) THEN
+                WHERE (contzones == contzones(i,ij+pzones(k)-1) .AND. contzones /= Ncontmin)
+                  contzones = Ncontmin
+                END WHERE
+                allzones(contzones(i,ij+pzones(k)-1)) = 0
+              END IF
+            END DO
+            Ncont = Ncontmin
+          END IF
+
+          IF (.NOT.assigned(i,j)) THEN
+            Nzones = Nzones + 1
+            Ncont = Nzones
+            allzones(Nzones) = 1
+          ELSE
+            Ncont = contzones(i,j)
+          END IF
+          contzones(i,j) = Ncont
+          contzones(ii:ei,j) = Ncont
+          contzones(i,ij:ej) = Ncont
+          notdone(i,j) = .FALSE.
+          assigned(i,j) = .TRUE.
+          assigned(ii:ei,j) = .TRUE.
+          assigned(i,ij:ej) = .TRUE.
+        END IF
+      END DO
+    END DO
+
+    ! Using allzones to provide continuous assigned values
+    Nzones = 0
+    DO k=1, dx*dy/4
+      IF (allzones(k) /= 0) THEN
+        Nzones = Nzones + 1
+        IF (allzones(k) /= Nzones) THEN
+          WHERE(contzones == allzones(k))
+            contzones = Nzones
+          END WHERE
+        END IF
+      END IF
+    END DO
+
+    RETURN
+
+  END SUBROUTINE continguos_homogene_zones
+
+  SUBROUTINE get_xyconlimits(d1, d2, matv, NOval, i, j, ix, ex, iy, ey)
+  ! Subroutine for getting the limits of contiguous values from a given point in a 2D matrix
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1, d2, i, j, NOval
+    INTEGER, DIMENSION(d1,d2), INTENT(in)                :: matv
+    INTEGER, INTENT(out)                                 :: ix, ex, iy, ey
+
+    ! Local
+    INTEGER                                              :: i1, j1
+    LOGICAL                                              :: found
+
+!!!!!!! Variables
+! d1, d2: Shape of input 2D values
+! i, j: grid point from which get the contiguous values
+! NOval: value for grid points without data
+! matv: 2D matrx values
+! ix, ex, iy, ey: limits from a grid point
+
+    fname = 'get_xyconlimits'
+
+    ix = -1
+    ex = -1
+    iy = -1
+    ey = -1
+
+    ! Before i
+    IF (i > 1) THEN
+      ix = i
+      found = .FALSE.
+      DO i1=i-1,1,-1
+        IF (matv(i1,j) == NOval) THEN
+          ix = i1 + 1
+          found = .TRUE.
+          EXIT
+        END IF
+      END DO
+      IF (.NOT.found) ix=1
+    ELSE
+      ix = 1
+    END IF
+    ! After i
+    IF (i < d1) THEN
+      ex = i
+      found = .FALSE.
+      DO i1=i+1,d1
+        IF (matv(i1,j) == NOval) THEN
+          ex = i1 - 1
+          found = .TRUE.
+          EXIT
+        END IF
+      END DO
+      IF (.NOT.found) ex=d1
+    ELSE
+      ex = d1
+    END IF
+
+    ! Before j
+    IF (j > 1) THEN
+      iy = j
+      found = .FALSE.
+      DO j1=j-1,1,-1
+        IF (matv(i,j1) == NOval) THEN
+          iy = j1 + 1
+          found = .TRUE.
+          EXIT
+        END IF
+      END DO
+      IF (.NOT.found) iy=1
+    ELSE
+      iy = 1
+    END IF
+    ! After j
+    IF (j < d2) THEN
+      ey = j
+      found = .FALSE.
+      DO j1=j+1,d2
+        IF (matv(i,j1) == NOval) THEN
+          ey = j1 - 1
+          found = .TRUE.
+          EXIT
+        END IF
+      END DO
+      IF (.NOT.found) ey=d2
+    ELSE
+      ey = d2
+    END IF
+
+    RETURN
+
+  END SUBROUTINE get_xyconlimits
+
 END MODULE module_generic