![(please configure the [header_logo] section in trac.ini)](/LMDZ_WRF/chrome/site/your_project_logo.png){kind=logo}
Changeset 1651 in lmdz_wrf for trunk/tools
- Timestamp:
- Sep 19, 2017, 6:10:56 PM (8 years ago)
- File:
-
- 1 edited
-
trunk/tools/module_scientific.f90 (modified) (19 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/tools/module_scientific.f90
r1650 r1651 37 37 CONTAINS 38 38 39 39 40 SUBROUTINE poly_overlap_tracks(dbg, dx, dy, dt, Nallpolys, allpolys, ctrpolys, areapolys, Nmaxpoly, & 40 41 Nmaxtracks, tracks, finaltracks) … … 46 47 INTEGER, INTENT(in) :: dx, dy, dt, Nmaxpoly, Nmaxtracks 47 48 INTEGER, DIMENSION(dt), INTENT(in) :: Nallpolys 48 INTEGER, DIMENSION(d t,dx,dy), INTENT(in) :: allpolys49 REAL(r_k), DIMENSION( dt,Nmaxpoly,2), INTENT(in) :: ctrpolys50 REAL(r_k), DIMENSION( dt,Nmaxpoly), INTENT(in) :: areapolys51 REAL(r_k), DIMENSION( dt,Nmaxtracks,Nmaxpoly,5), &49 INTEGER, DIMENSION(dx,dy,dt), INTENT(in) :: allpolys 50 REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in) :: ctrpolys 51 REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in) :: areapolys 52 REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,dt), & 52 53 INTENT(out) :: tracks 53 REAL(r_k), DIMENSION( dt,Nmaxtracks,4), INTENT(out) :: finaltracks54 REAL(r_k), DIMENSION(4,Nmaxtracks,dt), INTENT(out) :: finaltracks 54 55 55 56 ! Local 56 57 INTEGER :: i, j, ip, it, iip, itt 57 58 INTEGER :: ierr 58 INTEGER, DIMENSION(dt,Nmaxpoly) :: coincidencies, NOcoincidencies, & 59 MULTIcoincidencies 59 INTEGER, DIMENSION(Nmaxpoly,dt) :: coincidencies, NOcoincidencies 60 60 INTEGER, DIMENSION(:), ALLOCATABLE :: coins 61 61 INTEGER, DIMENSION(:,:), ALLOCATABLE :: coinsNpts 62 INTEGER, DIMENSION( dt,Nmaxpoly) :: polycoincidencies63 INTEGER, DIMENSION( dt,Nmaxpoly,Nmaxpoly) :: coincidenciesNpts62 INTEGER, DIMENSION(Nmaxpoly,dt) :: polycoincidencies 63 INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt) :: coincidenciesNpts 64 64 INTEGER :: Nmeet, Nsearch, Nindep 65 65 INTEGER, DIMENSION(dt) :: Nindeppolys 66 66 CHARACTER(len=5) :: NcoinS 67 INTEGER, DIMENSION( dt,Nmaxpoly,Nmaxpoly) :: polysIndep68 INTEGER, DIMENSION( dt,Nmaxpoly) :: NpolysIndep69 INTEGER, DIMENSION( dt,Nmaxpoly) :: SpolysIndep67 INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt) :: polysIndep 68 INTEGER, DIMENSION(Nmaxpoly,dt) :: NpolysIndep 69 INTEGER, DIMENSION(Nmaxpoly,dt) :: SpolysIndep 70 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 71 INTEGER :: Nprev, NNprev, Ntprev 72 LOGICAL :: Indeppolychained 77 73 INTEGER :: itrack, ictrack 78 74 INTEGER :: ixp, iyp, ttrack 79 INTEGER :: Nnew80 LOGICAL, DIMENSION(Nmaxpoly) :: chained81 75 INTEGER, DIMENSION(dt) :: Ntracks 82 76 INTEGER :: idtrack, maxtrack 83 INTEGER, DIMENSION(dt,Nmaxpoly) :: tracks_it84 77 85 78 !!!!!!! Variables … … 104 97 ! Polygons without a coincidence 105 98 NOcoincidencies = 0 106 ! Polygons with multiple coincidences107 MULTIcoincidencies = 0108 99 ! Number of independent polygons by time step 109 100 Nindeppolys = 0 … … 120 111 Nindeppolys(it) = Nmeet 121 112 DO i=1, Nmeet 122 SpolysIndep(i t,i) = i123 NpolysIndep( it,1:Nmeet) = 1124 polysIndep( it,i,1) = i113 SpolysIndep(i,it) = i 114 NpolysIndep(1:Nmeet,it) = 1 115 polysIndep(1,i,it) = i 125 116 END DO 126 117 … … 140 131 CALL ErrMsg(msg,fname,ierr) 141 132 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) = coins146 coincidenciesNpts( it+1,1:Nmeet,1:Nsearch) = coinsNpts133 CALL coincidence_all_polys(dbg, dx, dy, Nmeet, allpolys(:,:,it+1), ctrpolys(:,1:Nmeet,it+1), & 134 Nsearch, allpolys(:,:,it), ctrpolys(:,1:Nsearch,it), areapolys(1:Nsearch,it), coins, coinsNpts) 135 136 polycoincidencies(1:Nmeet,it+1) = coins 137 coincidenciesNpts(1:Nmeet,1:Nsearch,it+1) = coinsNpts 147 138 148 139 ! Counting the number of times that a polygon has a coincidency … … 158 149 IF (coins(i) == -1) THEN 159 150 Nindep = Nindep + 1 160 NOcoincidencies(i t+1,i) = 1161 SpolysIndep( it+1,Nindep) = -1162 NpolysIndep( it+1,Nindep) = NpolysIndep(it+1,Nindep) + 1163 polysIndep( it+1,Nindep,NpolysIndep(it+1,Nindep)) = i151 NOcoincidencies(i,it+1) = 1 152 SpolysIndep(Nindep,it+1) = -1 153 NpolysIndep(Nindep,it+1) = NpolysIndep(Nindep,it+1) + 1 154 polysIndep(Nindep,NpolysIndep(Nindep,it+1),it+1) = i 164 155 ELSE IF (coins(i) == -9) THEN 165 156 WRITE(NcoinS,'(I5)')coins(i) … … 170 161 DO ip=1, Nsearch 171 162 IF (coins(i) == ip) THEN 172 IF (coincidencies(i t+1,ip) == 0) THEN163 IF (coincidencies(ip,it+1) == 0) THEN 173 164 Nindep = Nindep + 1 174 SpolysIndep( it+1,Nindep) = ip165 SpolysIndep(Nindep,it+1) = ip 175 166 END IF 176 coincidencies(i t+1,ip) = coincidencies(it+1,ip) + 1167 coincidencies(ip,it+1) = coincidencies(ip,it+1) + 1 177 168 DO iindep=1, Nindep 178 IF (SpolysIndep(i t+1,iindep) == coins(i)) THEN179 NpolysIndep(i t+1,iindep) = NpolysIndep(it+1,iindep) + 1180 polysIndep(i t+1,iindep,NpolysIndep(it+1,iindep)) = i169 IF (SpolysIndep(iindep,it+1) == coins(i)) THEN 170 NpolysIndep(iindep,it+1) = NpolysIndep(iindep,it+1) + 1 171 polysIndep(iindep,NpolysIndep(iindep,it+1),it+1) = i 181 172 END IF 182 173 END DO … … 192 183 PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______' 193 184 DO i=1, Nindeppolys(it+1) 194 PRINT *,i, SpolysIndep(i t+1,i), NpolysIndep(it+1,i), ':', &195 polysIndep(i t+1,i,1:NpolysIndep(it+1,i))185 PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', & 186 polysIndep(i,1:NpolysIndep(i,it+1),it+1) 196 187 END DO 197 188 END IF … … 204 195 PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs' 205 196 DO ip=1, Nindeppolys(it) 206 PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(i t,ip), ':', &207 polysIndep(i t,ip,1:NpolysIndep(it,ip))197 PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':', & 198 polysIndep(ip,1:NpolysIndep(ip,it),it) 208 199 END DO 209 200 END DO … … 221 212 DO ip=1, Nindeppolys(1) 222 213 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 214 tracks(1,1,itrack,1) = itrack*1. 215 tracks(2,1,itrack,1) = SpolysIndep(ip,1) 216 tracks(3,1,itrack,1) = ctrpolys(1,ip,1) 217 tracks(4,1,itrack,1) = ctrpolys(2,ip,1) 218 tracks(5,1,itrack,1) = 1 229 219 Ntracks(1) = Ntracks(1) + 1 230 220 END DO … … 235 225 ! Indep polygons current time-step 236 226 current_poly: DO i=1, Nindeppolys(it) 237 IF (dbg) PRINT *,' curent poly:', i, 'Prev poly:', SpolysIndep(i t,i), ' N ass. polygons:', &238 NpolysIndep(i t,i), 'ass. poly:', polysIndep(it,i,1:NpolysIndep(it,i))227 IF (dbg) PRINT *,' curent poly:', i, 'Prev poly:', SpolysIndep(i,it), ' N ass. polygons:', & 228 NpolysIndep(i,it), 'ass. poly:', polysIndep(i,1:NpolysIndep(i,it),it) 239 229 Indeppolychained = .FALSE. 240 230 … … 242 232 ! Looping overall 243 233 it1_tracks: DO itt=1, Ntracks(it-1) 244 itrack = tracks( it-1,itt,1,1)234 itrack = tracks(1,1,itt,it-1) 245 235 ! 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)236 Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0) 237 IF (dbg) PRINT *,itt,' track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1) 248 238 249 239 ! Looking for coincidencies 250 240 DO iip=1, Ntprev 251 IF (tracks( it-1,itt,iip,2) == SpolysIndep(it,i)) THEN241 IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN 252 242 Indeppolychained = .TRUE. 253 IF (dbg) PRINT *,' coincidence found by polygon:', tracks( it-1,itt,iip,2)243 IF (dbg) PRINT *,' coincidence found by polygon:', tracks(2,iip,itt,it-1) 254 244 EXIT 255 245 END IF … … 259 249 ictrack = Ntracks(it) 260 250 ! Assigning all the IDs to the next step of the track 261 DO iip=1, NpolysIndep(i t,i)262 iiprev = polysIndep(i t,i,iip)263 tracks( it,ictrack,iip,1) = itrack264 tracks( it,ictrack,iip,2) = iiprev265 ixp = ctrpolys( it,iiprev,1)266 iyp = ctrpolys( it,iiprev,2)267 tracks( it,ictrack,iip,3) = ixp268 tracks( it,ictrack,iip,4) = iyp269 tracks( it,ictrack,iip,5) = it251 DO iip=1, NpolysIndep(i,it) 252 iiprev = polysIndep(i,iip,it) 253 tracks(1,iip,ictrack,it) = itrack 254 tracks(2,iip,ictrack,it) = iiprev 255 ixp = ctrpolys(1,iiprev,it) 256 iyp = ctrpolys(2,iiprev,it) 257 tracks(3,iip,ictrack,it) = ixp 258 tracks(4,iip,ictrack,it) = iyp 259 tracks(5,iip,ictrack,it) = it 270 260 END DO 271 tracks_it(it,ictrack) = itrack272 261 EXIT 273 262 END IF … … 279 268 ictrack = Ntracks(it) 280 269 ! ID of new track 281 maxtrack = INT(MAXVAL(tracks( :,:,:,1)*1.))270 maxtrack = INT(MAXVAL(tracks(1,:,:,:)*1.)) 282 271 IF (dbg) PRINT *,' New track!', maxtrack+1 283 272 284 273 ! Assigning all the IDs to the next step of the track 285 DO j=1, NpolysIndep(i t,i)286 iiprev = polysIndep(i t,i,j)287 tracks( it,ictrack,j,1) = maxtrack+1288 tracks( it,ictrack,j,2) = iiprev289 ixp = ctrpolys( it,iiprev,1)290 iyp = ctrpolys( it,iiprev,2)291 tracks( it,ictrack,j,3) = ixp292 tracks( it,ictrack,j,4) = iyp293 tracks( it,ictrack,j,5) = it274 DO j=1, NpolysIndep(i,it) 275 iiprev = polysIndep(i,j,it) 276 tracks(1,j,ictrack,it) = maxtrack+1 277 tracks(2,j,ictrack,it) = iiprev 278 ixp = ctrpolys(1,iiprev,it) 279 iyp = ctrpolys(2,iiprev,it) 280 tracks(3,j,ictrack,it) = ixp 281 tracks(4,j,ictrack,it) = iyp 282 tracks(5,j,ictrack,it) = it 294 283 END DO 295 tracks_it(it,Ntracks(it)) = maxtrack+1296 284 END IF 297 285 … … 300 288 IF (dbg) PRINT *,' At this time-step:', it, ' N trajectories:', Ntracks(it) 301 289 DO i=1, Ntracks(it) 302 Nprev = COUNT(INT(tracks( it,i,:,2)) /= 0)303 PRINT *,i,'tracks:', tracks( it,i,1:Nprev,2)290 Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0) 291 PRINT *,i,'tracks:', tracks(2,1:Nprev,i,it) 304 292 END DO 305 293 … … 310 298 311 299 finaltracks = 0. 312 maxtrack = MAXVAL(tracks( :,:,:,1))300 maxtrack = MAXVAL(tracks(1,:,:,:)) 313 301 314 302 DO it=1, dt 315 303 DO itt=1, Ntracks(it) 316 itrack = INT(tracks( it,itt,1,1))317 Nprev = COUNT(INT(tracks( it,itt,:,2)) /= 0)304 itrack = INT(tracks(1,1,itt,it)) 305 Nprev = COUNT(INT(tracks(2,:,itt,it)) /= 0) 318 306 PRINT *,'it:', it,'itrack:', itrack, 'Nprev:', Nprev 319 finaltracks( it,itrack,1) = itrack*1.320 finaltracks( it,itrack,2) = SUM(tracks(it,itt,:,3))/Nprev321 finaltracks( it,itrack,3) = SUM(tracks(it,itt,:,4))/Nprev322 finaltracks( it,itrack,4) = it*1.323 PRINT *,' finaltrack:', finaltracks( it,itrack,:)307 finaltracks(1,itrack,it) = itrack*1. 308 finaltracks(2,itrack,it) = SUM(tracks(3,:,itt,it))/Nprev 309 finaltracks(3,itrack,it) = SUM(tracks(4,:,itt,it))/Nprev 310 finaltracks(4,itrack,it) = it*1. 311 PRINT *,' finaltrack:', finaltracks(:,itrack,it) 324 312 END DO 325 313 END DO … … 339 327 INTEGER, INTENT(in) :: dx, dy, Nallpoly, Npoly 340 328 INTEGER, DIMENSION(dx,dy), INTENT(in) :: allpoly, polys 341 REAL(r_k), DIMENSION( Nallpoly,2), INTENT(in) :: icpolys342 REAL(r_k), DIMENSION( Npoly,2), INTENT(in) :: cpolys329 REAL(r_k), DIMENSION(2,Nallpoly), INTENT(in) :: icpolys 330 REAL(r_k), DIMENSION(2,Npoly), INTENT(in) :: cpolys 343 331 REAL(r_k), DIMENSION(Npoly), INTENT(in) :: apolys 344 332 INTEGER, DIMENSION(Nallpoly), INTENT(out) :: polycoins … … 383 371 DO j=1, Npoly 384 372 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 )373 dist = SQRT( (icpolys(1,ip)*1.-cpolys(1,j)*1.)**2 + (icpolys(2,ip)*1.-cpolys(2,j)*1.)**2 ) 386 374 IF ( dist > maxcoindist) THEN 387 375 maxcoindist = dist … … 481 469 482 470 END SUBROUTINE coincidence_poly 483 484 471 485 472 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.
