Changeset 2332 in lmdz_wrf

Timestamp:

Feb 12, 2019, 2:41:13 PM (6 years ago)

Author:

lfita

Message:

Adding use of 'polygons' at the end of `rangefaces' in order to homogenously join ranges

Location:

trunk/tools

Files:

• diag_tools.py (modified) (1 diff)
• module_ForDiagnostics.f90 (modified) (3 diffs)
• module_scientific.f90 (modified) (20 diffs)

trunk/tools/diag_tools.py

@@ -1480 +1480 @@
     """
     fname = 'Forcompute_range_faces'
-
-    print fname + 'Lluis:', dsfilt, dsnewrng, hvalleyrng
 
     vardims = dimns[:]

trunk/tools/module_ForDiagnostics.f90

@@ -9 +9 @@
   USE module_definitions
   USE module_generic
+  USE module_scientific
   USE module_ForDiagnosticsVars
 
@@ -978 +979 @@
     INTEGER, DIMENSION(2,d2)                             :: ranges2
     INTEGER, DIMENSION(d1,d2)                            :: iranges
+    LOGICAL, DIMENSION(d1,d2)                            :: Lranges
 
 !!!!!!! Variables
@@ -1053 +1055 @@
 
     ! Homogenizing indices of the ranges
-    IF (TRIM(face) == 'WE') THEN
-      CALL xzones_homogenization(d1, d2, iranges, ranges)
-    ELSE IF (TRIM(face) == 'SN') THEN
-      CALL yzones_homogenization(d1, d2, iranges, ranges)
-    END IF
+    Lranges = iranges /= 0
+    CALL polygons(.FALSE., d1, d2, Lranges, ranges, Nranges)
+
+!    IF (TRIM(face) == 'WE') THEN
+!      CALL xzones_homogenization(d1, d2, iranges, ranges)
+!    ELSE IF (TRIM(face) == 'SN') THEN
+!      CALL yzones_homogenization(d1, d2, iranges, ranges)
+!    END IF
 
     RETURN

trunk/tools/module_scientific.f90

@@ -2576 +2576 @@
   INTEGER                                                :: ierr
   INTEGER, DIMENSION(:,:), ALLOCATABLE                   :: borders
-  LOGICAL, DIMENSION(dx,dy)                              :: isborder, isbordery
+  LOGICAL, DIMENSION(dx,dy)                              :: isborder, isbordery, borderp
   INTEGER, DIMENSION(:,:,:), ALLOCATABLE                 :: paths
   INTEGER                                                :: Npath
@@ -2584 +2584 @@
   INTEGER, DIMENSION(:,:), ALLOCATABLE                   :: vertxs, points
   LOGICAL, DIMENSION(:), ALLOCATABLE                     :: isin
+  CHARACTER(len=1000)                                    :: boundsS
 
 !!!!!!! Variables
@@ -2596 +2597 @@
 
   ! The mathematical maximum woiuld be dx*dy/4, but let's be optimistic... (sorry Jero)
-  Nppt = dx*dy/10
+  Nppt = dx*dy/100
 
   IF (ALLOCATED(borders)) DEALLOCATE(borders)
@@ -2605 +2606 @@
   IF (ALLOCATED(paths)) DEALLOCATE(paths)
   ALLOCATE(paths(Nppt,Nppt,2), STAT=ierr)
-  msg = "Problems allocating matrix 'paths'"
+  boundsS = vectorI_S(3, (/Nppt, Nppt, 2/))
+  msg = "Problems allocating matrix 'paths' shape: " // TRIM(boundsS) // " try to reduce Nppt " //    &
+    "and recompile"
   CALL ErrMsg(msg, fname, ierr)
 
@@ -2633 +2636 @@
   END DO
 
-  CALL borders_matrixL(dx, dy, Nppt, boolmat, borders, isborder, isbordery)
+  CALL borders_matrixL(dbg, dx, dy, Nppt, boolmat, borders, isborder, isbordery)
   CALL paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)
 
@@ -2651 +2654 @@
     isin = .FALSE.
 
-    IF (dbg) PRINT *, '  path:', ip, ' N pts:', Nptpaths(ip)
+    IF (dbg) THEN
+      PRINT *, '  path:', ip, ' N pts:', Nptpaths(ip)
+      DO j=1, Nptpaths(ip)
+        PRINT *, '      ',j,':',paths(ip,j,:)
+      END DO
+    END IF
+
+    borderp = .FALSE.
+    DO j=1,Nptpaths(ip)
+      borderp(paths(ip,j,1),paths(ip,j,2)) = .TRUE.
+    END DO
 
     CALL path_properties(dx, dy, boolmat, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), xtrx, xtry,       &
@@ -2667 +2680 @@
     END IF
  
-    CALL gridpoints_InsidePolygon(dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), Nvertx,&
-      xtrx, xtry, vertxs, Npts, points, isin)
+    CALL gridpoints_InsidePolygon(dbg, dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:),   &
+      Nvertx, xtrx, xtry, vertxs, Npts, points, isin)
 
     ! Filling polygons
@@ -2689 +2702 @@
   CALL clean_polygons(dx, dy, boolmat, polys, Npoly, dbg)
 
-  DEALLOCATE (borders)  
-  DEALLOCATE (Nptpaths)
-  DEALLOCATE (paths)
-  DEALLOCATE (vertxs)
-  DEALLOCATE (points)
-  DEALLOCATE (isin)
+  IF (ALLOCATED(borders)) DEALLOCATE (borders)  
+  IF (ALLOCATED(Nptpaths)) DEALLOCATE (Nptpaths)
+  IF (ALLOCATED(paths)) DEALLOCATE (paths)
+  IF (ALLOCATED(vertxs)) DEALLOCATE (vertxs)
+  IF (ALLOCATED(points)) DEALLOCATE (points)
+  IF (ALLOCATED(isin)) DEALLOCATE (isin)
 
   RETURN
@@ -2904 +2917 @@
   END SUBROUTINE path_properties
 
-  SUBROUTINE gridpoints_InsidePolygon(dx, dy, isbrdr, Npath, path, Nvrtx, xpathxtrm, ypathxtrm,       &
+  SUBROUTINE gridpoints_InsidePolygon(dbg, dx, dy, isbrdr, Npath, path, Nvrtx, xpathxtrm, ypathxtrm,  &
     vrtxs, Npts, pts, inside)
 ! Subroutine to determine if a series of grid points are inside a polygon following ray casting algorithm
@@ -2912 +2925 @@
 
   INTEGER, INTENT(in)                                    :: dx,dy,Npath,Nvrtx,Npts
+  LOGICAL, INTENT(in)                                    :: dbg
   LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: isbrdr
   INTEGER, DIMENSION(Npath,2), INTENT(in)                :: path
@@ -2947 +2961 @@
   halo_brdr = .FALSE.
 
+  IF (dbg) PRINT *,'Border _______'
   DO i=1,dx
     halo_brdr(i+1,2:dy+1) = isbrdr(i,:)
+    IF (dbg) PRINT *,isbrdr(i,:)
   END DO
 
@@ -2979 +2995 @@
         IF (halo_brdr(ix+1,j+1) .AND. (ispath /= -1) .AND. (halo_brdr(ix+1,j+1) .EQV. halo_brdr(ix+1,j+2))) THEN
           Nbrbrdr = Nbrbrdr + 1
+          IF (dbg) PRINT *,'    ',Nbrbrdr,' Consec brdrs:', halo_brdr(ix+1,j+1), halo_brdr(ix+1,j+2), &
+             '(', ix,j,';', ix,j+1,')', isbrdr(ix,j), isbrdr(ix,j+1) 
         ELSE
           ! Will remove that consecutive borders above 2
           IF (Nbrbrdr /= 0) THEN
+            IF (dbg) PRINT *, ix,',',iy,';', Nintersecs, '  amount of consecutive borders:', Nbrbrdr, &
+              ' removing:', MAX(Nbrbrdr-1, 0)
             Nintersecs = Nintersecs - MAX(Nbrbrdr-1, 0)
             Nbrbrdr = 0
@@ -2988 +3008 @@
       END DO
       IF (MOD(Nintersecs,2) /= 0) inside(ip) = .TRUE.
+      IF (dbg) PRINT *,ip,'    point:', ix, iy, 'isbrdr:', isbrdr(ix,1:iy-1), 'y-ray:', halo_brdr(ix+1,1:iy), 'inside:', inside(ip)
     END IF
 
@@ -3065 +3086 @@
 END SUBROUTINE look_clockwise_borders
 
-SUBROUTINE borders_matrixL(dx,dy,dxy,Lmat,brdrs,isbrdr,isbrdry)
+SUBROUTINE borders_matrixL(dbg,dx,dy,dxy,Lmat,brdrs,isbrdr,isbrdry)
 ! Subroutine to provide the borders of a logical array (interested in .TRUE.)
 
@@ -3071 +3092 @@
 
   INTEGER, INTENT(in)                                    :: dx,dy,dxy
+  LOGICAL, INTENT(in)                                    :: dbg
   LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: Lmat
   INTEGER, DIMENSION(dxy,2), INTENT(out)                 :: brdrs
@@ -3185 +3207 @@
   END DO
 
+  IF (dbg) THEN
+    PRINT *,' BORDERS _______ x y'
+    DO i=1,dx
+      PRINT *,isbrdr(i,:), '       ', isbrdry(i,:)
+    END DO
+  END IF
+
   RETURN
 
@@ -3222 +3251 @@
   fname = 'paths_border'
 
+  IF (dbg) PRINT *, TRIM(fname) // ' ...'
+
   ! Sarting matrix
   paths = -1
@@ -3238 +3269 @@
     
   IF (dbg) THEN
+    PRINT *,'  isborder ______'
+    DO i=1,dx
+      PRINT *,isborder(i,:)
+    END DO
+
     PRINT *,'    borders _______'
     DO i=1,Nbrdr
@@ -3278 +3314 @@
         IF (dbg) PRINT *,ipath,'iip jip:', iip, ijp
         found = .FALSE.
-        CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE.,iif,jjf,iff)
+        CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,dbg,iif,jjf,iff)
         IF (iif /= -1) THEN
           ip=ip+1
@@ -3310 +3346 @@
               iip = paths(ipath,i,1)
               ijp = paths(ipath,i,2)
-              CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE., iif,   &
+              CALL look_clockwise_borders(dx,dy,Nppt,borders, gotbrdr, isborder,iip, ijp, dbg, iif,   &
                 jjf,iff)
               IF (iif /= -1 .AND. iff /= -1) THEN