Changeset 1654 in lmdz_wrf

Timestamp:

Sep 25, 2017, 11:28:03 PM (7 years ago)

Author:

lfita

Message:

Fixing:

• calcualtion of mean weighted-center of the polygon at `polygon_properties'
• Adding compilation at single resolution

Location:

trunk/tools

Files:

• interpolate.f90 (modified) (1 diff)
• module_definitions.f90 (modified) (1 diff)
• module_scientific.f90 (modified) (15 diffs)

trunk/tools/interpolate.f90

@@ -74 +74 @@
   difffraclonlat = SQRT((fraclon-lonv)**2. + (fraclat-latv)**2.)
   mindifffracLl = MINVAL(difffraclonlat)
+!  ilonlatfrac = index2DArrayR_K(difffraclonlat, Nperx, Npery, mindifffracLl)
   ilonlatfrac = index2DArrayR(difffraclonlat, Nperx, Npery, mindifffracLl)
 

trunk/tools/module_definitions.f90

@@ -1 +1 @@
 MODULE module_definitions
 
-  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+!  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
+  INTEGER, PARAMETER                                     :: r_k = KIND(1.)
   REAL(r_k), PARAMETER                                   :: zeroRK = 0.
   REAL(r_k), PARAMETER                                   :: oneRK = 1.

trunk/tools/module_scientific.f90

@@ -43 +43 @@
   CONTAINS
 
-
   SUBROUTINE poly_overlap_tracks_area(dbg, dx, dy, dt, minarea, inNallpolys, allpolys, ctrpolys,      &
     areapolys, Nmaxpoly, Nmaxtracks, tracks, finaltracks)
@@ -65 +64 @@
     INTEGER                                              :: i, j, ip, it, iip, itt
     INTEGER                                              :: ierr
-    REAL(r_k), DIMENSION(2,Nmaxpoly,dt)                  :: Actrpolys
-    REAL(r_k), DIMENSION(Nmaxpoly,dt)                    :: Aareapolys
+    REAL(r_k), DIMENSION(Nmaxpoly)                       :: Aprevpolys, Acurrpolys
+    REAL(r_k), DIMENSION(2,Nmaxpoly)                     :: Cprevpolys, Ccurrpolys
     INTEGER, DIMENSION(dt)                               :: Nallpolys
     INTEGER, DIMENSION(dx,dy)                            :: meetpolys, searchpolys
@@ -73 +72 @@
     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
@@ -105 +102 @@
     IF (dbg) PRINT *,TRIM(fname)
 
-    polycoincidencies = fillvalI
-    coincidenciesNpts = fillvalI
     ! Number of independent polygons by time step
     Nindeppolys = 0
@@ -131 +126 @@
         SpolysIndep(ip,it) = i
         currID(ip) = i
+        Acurrpolys(ip) = areapolys(i,it)
+        Ccurrpolys(1,ip) = ctrpolys(1,i,it)
+        Ccurrpolys(2,ip) = ctrpolys(2,i,it)
         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)
@@ -167 +162 @@
       prevID = 0
       prevID = currID
+      Aprevpolys = Acurrpolys
+      Cprevpolys = Ccurrpolys
       ip = 0
 
@@ -173 +170 @@
           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)
+          Acurrpolys(ip) = areapolys(i,it+1)
+          Acurrpolys(ip) = areapolys(i,it+1)
+          Ccurrpolys(1,ip) = ctrpolys(1,i,it+1)
+          Ccurrpolys(2,ip) = ctrpolys(2,i,it+1)
         ELSE
           WHERE (meetpolys == i)
@@ -199 +197 @@
       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,     &
+      CALL coincidence_all_polys_area(dbg, dx,dy, Nmeet, currID, meetpolys, Acurrpolys(1:Nmeet),      &
+        Nsearch, prevID, searchpolys, Cprevpolys(:,1:Nsearch), Aprevpolys(1:Nsearch), 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
@@ -234 +229 @@
             DO j=1, NNprev
               IF (coins(i) == polysIndep(ip,j,it)) THEN
-                IF (coincidencies(ip) == 0) THEN
+                ! Which index corresponds to this coincidence?
+                iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
+                IF (iindep == -1) 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))
+                iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
                 IF (iindep < 0) THEN
+                  PRINT *,'    Looking for:', coins(i)
+                  PRINT *,'    Within:', SpolysIndep(1:Nindep,it+1)
+                  PRINT *,'    Might content:', polysIndep(ip,1:NNprev,it)
+                  PRINT *,'    From an initial list:', coins(1:Nmeet)
                   msg = 'Wrong index! There must be an index here'
                   CALL ErrMsg(msg,fname,iindep)
@@ -384 +384 @@
         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
+
+    DEALLOCATE(coins)
+    DEALLOCATE(coinsNpts)
 
     RETURN
@@ -1222 +1219 @@
   meanvnctr = 0.
   meanvxctr = 0.
+  xcorr = 0.
+  ycorr = 0.
 
   nval = fillval64
@@ -1240 +1239 @@
           CASE ('lon/lat')
             ! Area as fixed yres and sinus-lat varying for xres
-            xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
+!            IF (KIND(xcorr(i,j)) == KIND(1.d0)) THEN
+!              xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
+!            ELSE
+              xcorr(i,j) = ABS(SIN(lon(i,j)*DegRad))
+!            END IF
             ycorr(i,j) = 1.
           CASE ('nadir-sat')
             ! Area from nadir resolution and degrading as we get far from satellite's nadir
             ! GOES-E: 0 N, 75 W
-            xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
+!            IF (KIND(xcorr(i,j)) == KIND(1.d0)) THEN
+!              xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
+!            ELSE
+              xcorr(i,j) = ABS(SIN(lon(i,j)*DegRad))
+!            END IF
             ycorr(i,j) = 1.
         END SELECT
@@ -1267 +1274 @@
 
   IF (dbg) THEN
-    PRINT *,'  grdi_coord lon lat value _______ '
+    PRINT *,'  grid_coord lon lat value _______ '
     DO ip=1, Npolypts
       PRINT *, polypts(ip,:), ';', lon(polypts(ip,1),polypts(ip,2)), lat(polypts(ip,1),polypts(ip,2)),&
@@ -1285 +1292 @@
     PRINT *,'  mean grid center: ', meanctr, ' weighted mean center: ', meanwctr
   END IF
-  
+
   ! Statistics of the values within the polygon
   CALL StatsR_K(Npolypts, polyvals, nval, xval, mval, m2val, stdval)
@@ -1296 +1303 @@
 
   ! Mean center weighted to minimum and maximum values
-  distvn = DABS(polyvals - nval)
-  distvx = DABS(polyvals - xval)
+!  IF (KIND(polyvals(1)) == KIND(1.d0)) THEN
+!    distvn = DABS(polyvals - nval)
+!    distvx = DABS(polyvals - xval)
+!  ELSE
+    distvn = ABS(polyvals - nval)
+    distvx = ABS(polyvals - xval)
+!  END IF
 
   meanvnctr = 0.