![(please configure the [header_logo] section in trac.ini)](/LMDZ_WRF/chrome/site/your_project_logo.png){kind=logo}
Changeset 1650 in lmdz_wrf
- Timestamp:
- Sep 18, 2017, 10:27:19 PM (7 years ago)
- File:
-
- 1 edited
-
trunk/tools/module_scientific.f90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/tools/module_scientific.f90
r1622 r1650 5 5 ! all_polygons_properties: Subroutine to determine the properties of all polygons in a 2D field: 6 6 ! borders_matrixL: Subroutine to provide the borders of a logical array (interested in .TRUE.) 7 ! coincidence_all_polys: Subtourine to determine which is the coincident polygon when a boolean polygon is provided 8 ! to a map of integer polygons 9 ! coincidence_poly: Subtourine to determine which is the coincident polygon when a boolean polygon is provided 10 ! to a map of integer polygons 7 11 ! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole 8 12 ! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End. … … 16 20 ! polygons: Subroutine to search the polygons of a border field. FORTRAN based. 1st = 1! 17 21 ! polygons_t: Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1! 18 ! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k) 22 ! poly_overlap_tracks: Subroutine to determine tracks of a series of consecutive 2D field with polygons using 23 ! maximum overlaping/coincidence! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k) 19 24 ! quantilesR_K: Subroutine to provide the quantiles of a given set of values of type real 'r_k' 20 25 ! rand_sample: Subroutine to randomly sample a range of indices … … 31 36 32 37 CONTAINS 38 39 SUBROUTINE poly_overlap_tracks(dbg, dx, dy, dt, Nallpolys, allpolys, ctrpolys, areapolys, Nmaxpoly, & 40 Nmaxtracks, tracks, finaltracks) 41 ! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum overlaping/coincidence 42 43 IMPLICIT NONE 44 45 LOGICAL, INTENT(in) :: dbg 46 INTEGER, INTENT(in) :: dx, dy, dt, Nmaxpoly, Nmaxtracks 47 INTEGER, DIMENSION(dt), INTENT(in) :: Nallpolys 48 INTEGER, DIMENSION(dt,dx,dy), INTENT(in) :: allpolys 49 REAL(r_k), DIMENSION(dt,Nmaxpoly,2), INTENT(in) :: ctrpolys 50 REAL(r_k), DIMENSION(dt,Nmaxpoly), INTENT(in) :: areapolys 51 REAL(r_k), DIMENSION(dt,Nmaxtracks,Nmaxpoly,5), & 52 INTENT(out) :: tracks 53 REAL(r_k), DIMENSION(dt,Nmaxtracks,4), INTENT(out) :: finaltracks 54 55 ! Local 56 INTEGER :: i, j, ip, it, iip, itt 57 INTEGER :: ierr 58 INTEGER, DIMENSION(dt,Nmaxpoly) :: coincidencies, NOcoincidencies, & 59 MULTIcoincidencies 60 INTEGER, DIMENSION(:), ALLOCATABLE :: coins 61 INTEGER, DIMENSION(:,:), ALLOCATABLE :: coinsNpts 62 INTEGER, DIMENSION(dt,Nmaxpoly) :: polycoincidencies 63 INTEGER, DIMENSION(dt,Nmaxpoly,Nmaxpoly) :: coincidenciesNpts 64 INTEGER :: Nmeet, Nsearch, Nindep 65 INTEGER, DIMENSION(dt) :: Nindeppolys 66 CHARACTER(len=5) :: NcoinS 67 INTEGER, DIMENSION(dt,Nmaxpoly,Nmaxpoly) :: polysIndep 68 INTEGER, DIMENSION(dt,Nmaxpoly) :: NpolysIndep 69 INTEGER, DIMENSION(dt,Nmaxpoly) :: SpolysIndep 70 INTEGER :: iindep, iiprev 71 INTEGER, DIMENSION(dt,Nmaxpoly,Nmaxpoly) :: polysIndepchains 72 INTEGER :: Nprev, polyprev, Nnotchained, Nprevchained 73 INTEGER :: NNprev, Ntprev 74 LOGICAL :: Indeppolychained, prevchained 75 INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly) :: polysIndepNotchained 76 INTEGER, DIMENSION(2) :: trackperiod 77 INTEGER :: itrack, ictrack 78 INTEGER :: ixp, iyp, ttrack 79 INTEGER :: Nnew 80 LOGICAL, DIMENSION(Nmaxpoly) :: chained 81 INTEGER, DIMENSION(dt) :: Ntracks 82 INTEGER :: idtrack, maxtrack 83 INTEGER, DIMENSION(dt,Nmaxpoly) :: tracks_it 84 85 !!!!!!! Variables 86 ! dx,dy,dt: space/time dimensions 87 ! Nallpolys: Vector with the number of polygons at each time-step 88 ! allpolys: Series of 2D field with the polygons 89 ! ctrpolys: center of the polygons 90 ! areapolys: area of the polygons 91 ! Nmaxpoly: Maximum possible number of polygons 92 ! Nmaxtracks: maximum number of tracks 93 ! tracks: series of consecutive polygons 94 ! trackperiod: period of the track in time-steps 95 96 fname = 'poly_overlap_tracks' 97 98 IF (dbg) PRINT *,TRIM(fname) 99 100 polycoincidencies = fillvalI 101 coincidenciesNpts = fillvalI 102 ! Number of times that a polygon has a coincidence 103 coincidencies = 0 104 ! Polygons without a coincidence 105 NOcoincidencies = 0 106 ! Polygons with multiple coincidences 107 MULTIcoincidencies = 0 108 ! Number of independent polygons by time step 109 Nindeppolys = 0 110 ! Number of polygons attached to each independent polygons by time step 111 NpolysIndep = 0 112 ! ID of searching polygon attached to each independent polygons by time step 113 SpolysIndep = 0 114 ! ID of polygons attached to each independent polygons by time step 115 polysIndep = 0 116 117 ! First time-step all are independent polygons 118 it = 1 119 Nmeet = Nallpolys(it) 120 Nindeppolys(it) = Nmeet 121 DO i=1, Nmeet 122 SpolysIndep(it,i) = i 123 NpolysIndep(it,1:Nmeet) = 1 124 polysIndep(it,i,1) = i 125 END DO 126 127 ! Looking for the coincidencies at each time step 128 DO it=1, dt-1 129 Nmeet = Nallpolys(it+1) 130 Nsearch = Nallpolys(it) 131 132 IF (ALLOCATED(coins)) DEALLOCATE(coins) 133 ALLOCATE(coins(Nmeet), STAT=ierr) 134 msg="Problems allocating 'coins'" 135 CALL ErrMsg(msg,fname,ierr) 136 137 IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts) 138 ALLOCATE(coinsNpts(Nmeet, Nsearch), STAT=ierr) 139 msg="Problems allocating 'coinsNpts'" 140 CALL ErrMsg(msg,fname,ierr) 141 142 CALL coincidence_all_polys(dbg, dx, dy, Nmeet, allpolys(it+1,:,:), ctrpolys(it+1,1:Nmeet,:), & 143 Nsearch, allpolys(it,:,:), ctrpolys(it,1:Nsearch,:), areapolys(it,1:Nsearch), coins, coinsNpts) 144 145 polycoincidencies(it+1,1:Nmeet) = coins 146 coincidenciesNpts(it+1,1:Nmeet,1:Nsearch) = coinsNpts 147 148 ! Counting the number of times that a polygon has a coincidency 149 IF (dbg) THEN 150 PRINT *,' Coincidencies for the given time-step:', it+1,' _______' 151 DO i=1, Nmeet 152 PRINT *,coins(i),' N search pts:', coinsNpts(i,:) 153 END DO 154 END IF 155 156 Nindep = 0 157 DO i=1, Nmeet 158 IF (coins(i) == -1) THEN 159 Nindep = Nindep + 1 160 NOcoincidencies(it+1,i) = 1 161 SpolysIndep(it+1,Nindep) = -1 162 NpolysIndep(it+1,Nindep) = NpolysIndep(it+1,Nindep) + 1 163 polysIndep(it+1,Nindep,NpolysIndep(it+1,Nindep)) = i 164 ELSE IF (coins(i) == -9) THEN 165 WRITE(NcoinS,'(I5)')coins(i) 166 msg="coins= "//TRIM(NcoinS)//" This is an error. One should have always only one " // & 167 "coincidence of polygon" 168 CALL ErrMsg(msg, fname, -1) 169 ELSE 170 DO ip=1, Nsearch 171 IF (coins(i) == ip) THEN 172 IF (coincidencies(it+1,ip) == 0) THEN 173 Nindep = Nindep + 1 174 SpolysIndep(it+1,Nindep) = ip 175 END IF 176 coincidencies(it+1,ip) = coincidencies(it+1,ip) + 1 177 DO iindep=1, Nindep 178 IF (SpolysIndep(it+1,iindep) == coins(i)) THEN 179 NpolysIndep(it+1,iindep) = NpolysIndep(it+1,iindep) + 1 180 polysIndep(it+1,iindep,NpolysIndep(it+1,iindep)) = i 181 END IF 182 END DO 183 END IF 184 END DO 185 END IF 186 END DO 187 Nindeppolys(it+1) = Nindep 188 189 IF (dbg) THEN 190 PRINT *,' time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',Nsearch 191 PRINT *,' number of independent polygons:', Nindeppolys(it+1) 192 PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______' 193 DO i=1, Nindeppolys(it+1) 194 PRINT *,i, SpolysIndep(it+1,i), NpolysIndep(it+1,i), ':', & 195 polysIndep(it+1,i,1:NpolysIndep(it+1,i)) 196 END DO 197 END IF 198 END DO 199 200 IF (dbg) THEN 201 PRINT *, 'Coincidencies to connect _______' 202 DO it=1, dt 203 PRINT *,' it:', it, ' Nindep:', Nindeppolys(it) 204 PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs' 205 DO ip=1, Nindeppolys(it) 206 PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(it,ip), ':', & 207 polysIndep(it,ip,1:NpolysIndep(it,ip)) 208 END DO 209 END DO 210 211 END IF 212 213 ! Trajectories 214 ! It should be done following the number of 'independent' polygons 215 ! One would concatenate that independent polygons which share IDs from one step to another 216 217 ! First time-step. Take all polygons 218 itrack = 0 219 tracks = 0. 220 Ntracks = 0 221 DO ip=1, Nindeppolys(1) 222 itrack = itrack + 1 223 tracks(1,itrack,1,1) = itrack*1. 224 tracks(1,itrack,1,2) = SpolysIndep(1,ip) 225 tracks(1,itrack,1,3) = ctrpolys(1,ip,1) 226 tracks(1,itrack,1,4) = ctrpolys(1,ip,2) 227 tracks(1,itrack,1,5) = 1 228 tracks_it(1,itrack) = itrack 229 Ntracks(1) = Ntracks(1) + 1 230 END DO 231 232 ! Looping allover already assigned tracks 233 timesteps: DO it=2, dt 234 IF (dbg) PRINT *,'timestep:', it, 'N indep polys:', Nindeppolys(it) 235 ! Indep polygons current time-step 236 current_poly: DO i=1, Nindeppolys(it) 237 IF (dbg) PRINT *,' curent poly:', i, 'Prev poly:', SpolysIndep(it,i), ' N ass. polygons:', & 238 NpolysIndep(it,i), 'ass. poly:', polysIndep(it,i,1:NpolysIndep(it,i)) 239 Indeppolychained = .FALSE. 240 241 ! Number of tracks previous time-step 242 ! Looping overall 243 it1_tracks: DO itt=1, Ntracks(it-1) 244 itrack = tracks(it-1,itt,1,1) 245 ! Number polygons ID assigned 246 Ntprev = COUNT(tracks(it-1,itt,:,2) /= 0) 247 IF (dbg) PRINT *,itt,' track:', itrack, 'assigned:', tracks(it-1,itt,1:Ntprev,2) 248 249 ! Looking for coincidencies 250 DO iip=1, Ntprev 251 IF (tracks(it-1,itt,iip,2) == SpolysIndep(it,i)) THEN 252 Indeppolychained = .TRUE. 253 IF (dbg) PRINT *,' coincidence found by polygon:', tracks(it-1,itt,iip,2) 254 EXIT 255 END IF 256 END DO 257 IF (Indeppolychained) THEN 258 Ntracks(it) = Ntracks(it) + 1 259 ictrack = Ntracks(it) 260 ! Assigning all the IDs to the next step of the track 261 DO iip=1, NpolysIndep(it,i) 262 iiprev = polysIndep(it,i,iip) 263 tracks(it,ictrack,iip,1) = itrack 264 tracks(it,ictrack,iip,2) = iiprev 265 ixp = ctrpolys(it,iiprev,1) 266 iyp = ctrpolys(it,iiprev,2) 267 tracks(it,ictrack,iip,3) = ixp 268 tracks(it,ictrack,iip,4) = iyp 269 tracks(it,ictrack,iip,5) = it 270 END DO 271 tracks_it(it,ictrack) = itrack 272 EXIT 273 END IF 274 END DO it1_tracks 275 276 ! Creation of a new track 277 IF (.NOT.Indeppolychained) THEN 278 Ntracks(it) = Ntracks(it) + 1 279 ictrack = Ntracks(it) 280 ! ID of new track 281 maxtrack = INT(MAXVAL(tracks(:,:,:,1)*1.)) 282 IF (dbg) PRINT *,' New track!', maxtrack+1 283 284 ! Assigning all the IDs to the next step of the track 285 DO j=1, NpolysIndep(it,i) 286 iiprev = polysIndep(it,i,j) 287 tracks(it,ictrack,j,1) = maxtrack+1 288 tracks(it,ictrack,j,2) = iiprev 289 ixp = ctrpolys(it,iiprev,1) 290 iyp = ctrpolys(it,iiprev,2) 291 tracks(it,ictrack,j,3) = ixp 292 tracks(it,ictrack,j,4) = iyp 293 tracks(it,ictrack,j,5) = it 294 END DO 295 tracks_it(it,Ntracks(it)) = maxtrack+1 296 END IF 297 298 END DO current_poly 299 300 IF (dbg) PRINT *,' At this time-step:', it, ' N trajectories:', Ntracks(it) 301 DO i=1, Ntracks(it) 302 Nprev = COUNT(INT(tracks(it,i,:,2)) /= 0) 303 PRINT *,i,'tracks:', tracks(it,i,1:Nprev,2) 304 END DO 305 306 END DO timesteps 307 308 ! Summarizing trajectories 309 ! When multiple polygons are available, the mean of their central positions determines the position 310 311 finaltracks = 0. 312 maxtrack = MAXVAL(tracks(:,:,:,1)) 313 314 DO it=1, dt 315 DO itt=1, Ntracks(it) 316 itrack = INT(tracks(it,itt,1,1)) 317 Nprev = COUNT(INT(tracks(it,itt,:,2)) /= 0) 318 PRINT *,'it:', it,'itrack:', itrack, 'Nprev:', Nprev 319 finaltracks(it,itrack,1) = itrack*1. 320 finaltracks(it,itrack,2) = SUM(tracks(it,itt,:,3))/Nprev 321 finaltracks(it,itrack,3) = SUM(tracks(it,itt,:,4))/Nprev 322 finaltracks(it,itrack,4) = it*1. 323 PRINT *,' finaltrack:', finaltracks(it,itrack,:) 324 END DO 325 END DO 326 327 RETURN 328 329 END SUBROUTINE poly_overlap_tracks 330 331 SUBROUTINE coincidence_all_polys(dbg, dx, dy, Nallpoly, allpoly, icpolys, Npoly, polys, cpolys, & 332 apolys, polycoins, coinNptss) 333 ! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons 334 ! In case of multiple coincidencies, the closest and then the biggest is taken 335 336 IMPLICIT NONE 337 338 LOGICAL, INTENT(in) :: dbg 339 INTEGER, INTENT(in) :: dx, dy, Nallpoly, Npoly 340 INTEGER, DIMENSION(dx,dy), INTENT(in) :: allpoly, polys 341 REAL(r_k), DIMENSION(Nallpoly,2), INTENT(in) :: icpolys 342 REAL(r_k), DIMENSION(Npoly,2), INTENT(in) :: cpolys 343 REAL(r_k), DIMENSION(Npoly), INTENT(in) :: apolys 344 INTEGER, DIMENSION(Nallpoly), INTENT(out) :: polycoins 345 INTEGER, DIMENSION(Nallpoly,Npoly), INTENT(out) :: coinNptss 346 347 ! Local 348 INTEGER :: i, j, ip 349 INTEGER :: maxcorr 350 INTEGER :: Nmaxcorr 351 LOGICAL, DIMENSION(dx,dy) :: boolpoly 352 INTEGER :: maxcoin 353 REAL :: dist, maxcoindist, maxcoinarea 354 355 !!!!!!! Variables 356 ! dx,dy: dimension of the space 357 ! Nallpoly: Number of polygons to find coincidence 358 ! allpoly: space with the polygons to meet 359 ! icpolys: center of the polygons to look for the coincidence 360 ! Npoly: number of polygons on the 2D space 361 ! polys: 2D field of polygons identified by their integer number (0 for no polygon) 362 ! cpolys: center of the polygons 363 ! apolys: area of the polygons 364 ! polycoins: coincident polyogn 365 ! -1: no-coincidence 366 ! 1 < Npoly: single coincidence with a given polygon 367 ! -9: coincidence with more than one polygon 368 ! coinNptss: number of points coincident with each polygon 369 370 fname = 'coincidence_all_polys' 371 IF (dbg) PRINT *,TRIM(fname) 372 373 DO ip=1, Nallpoly 374 boolpoly = allpoly == ip 375 CALL coincidence_poly(dbg, dx, dy, boolpoly, Npoly, polys, polycoins(ip), coinNptss(ip,:)) 376 IF (dbg) PRINT *,' polygon', ip, ' coincidence with:', polycoins(ip) 377 378 ! Coincidence with more than one polygon 379 IF (polycoins(ip) == -9) THEN 380 maxcoindist = -10. 381 maxcoinarea = -10. 382 maxcoin = MAXVAL(coinNptss(ip,:)) 383 DO j=1, Npoly 384 IF (coinNptss(ip,j) == maxcoin) THEN 385 dist = SQRT( (icpolys(ip,1)*1.-cpolys(j,1)*1.)**2 + (icpolys(ip,2)*1.-cpolys(j,2)*1.)**2 ) 386 IF ( dist > maxcoindist) THEN 387 maxcoindist = dist 388 maxcoinarea = apolys(j) 389 polycoins(ip) = j 390 ELSE IF ( dist == maxcoindist) THEN 391 IF (apolys(j) > maxcoinarea) THEN 392 polycoins(ip) = j 393 maxcoinarea = apolys(j) 394 END IF 395 END IF 396 END IF 397 END DO 398 END IF 399 END DO 400 401 RETURN 402 403 END SUBROUTINE coincidence_all_polys 404 405 SUBROUTINE coincidence_poly(dbg, dx, dy, poly, Npoly, polys, polycoin, coinNpts) 406 ! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons 407 408 IMPLICIT NONE 409 410 LOGICAL, INTENT(in) :: dbg 411 INTEGER, INTENT(in) :: dx, dy, Npoly 412 LOGICAL, DIMENSION(dx,dy), INTENT(in) :: poly 413 INTEGER, DIMENSION(dx,dy), INTENT(in) :: polys 414 INTEGER, INTENT(out) :: polycoin 415 INTEGER, DIMENSION(Npoly), INTENT(out) :: coinNpts 416 417 ! Local 418 INTEGER :: i, j, ip 419 INTEGER :: maxcorr 420 INTEGER :: Nmaxcorr 421 422 !!!!!!! Variables 423 ! dx,dy: dimension of the space 424 ! poly: bolean polygon to meet 425 ! Npoly: number of polygons on the 2D space 426 ! polys: 2D field of polygons identified by their integer number (0 for no polygon) 427 ! polycoin: coincident polyogn 428 ! -1: no-coincidence 429 ! 1 < Npoly: single coincidence with a given polygon 430 ! -9: coincidence with more than one polygon 431 ! coinNpts: number of points coincident with each polygon 432 433 fname = 'coincidence_poly' 434 IF (dbg) PRINT *,TRIM(fname) 435 436 IF (dbg) THEN 437 PRINT *,' Boolean polygon to search coincidences ...' 438 DO i=1,dx 439 PRINT *,poly(i,:) 440 END DO 441 442 PRINT *,' 2D polygons space ...' 443 DO i=1,dx 444 PRINT '(1000(I3,1x))',polys(i,:) 445 END DO 446 END IF 447 448 ! Looking for coincient points for the polygon 449 coinNpts = 0 450 DO i=1,dx 451 DO j=1,dy 452 IF (poly(i,j) .AND. polys(i,j) .NE. 0) coinNpts(polys(i,j)) = coinNpts(polys(i,j)) + 1 453 END DO 454 END DO 455 456 maxcorr = 0 457 maxcorr = INT(MAXVAL(coinNpts*1.)) 458 459 IF (dbg) PRINT *,' Maximum coincidence:', maxcorr 460 IF (maxcorr == 0) THEN 461 polycoin = -1 462 ELSE 463 Nmaxcorr = 0 464 DO ip=1, Npoly 465 IF (coinNpts(ip) == maxcorr) THEN 466 Nmaxcorr=Nmaxcorr+1 467 polycoin = ip 468 END IF 469 END DO 470 IF (Nmaxcorr > 1) polycoin = -9 471 END IF 472 473 IF (dbg) THEN 474 PRINT *,' Coincidences for each polygon _______' 475 DO ip=1, Npoly 476 PRINT *,' ', ip,': ', coinNpts(ip) 477 END DO 478 END IF 479 480 RETURN 481 482 END SUBROUTINE coincidence_poly 483 33 484 34 485 SUBROUTINE all_polygons_properties(dbg, dx, dy, Npoly, polys, lon, lat, values, xres, yres, projN, &
Note: See TracChangeset
for help on using the changeset viewer.
