![(please configure the [header_logo] section in trac.ini)](/LMDZ_WRF/chrome/site/your_project_logo.png){kind=logo}
Changeset 1660 in lmdz_wrf
- Timestamp:
- Oct 3, 2017, 8:09:21 PM (8 years ago)
- Location:
- trunk/tools
- Files:
-
- 5 edited
-
Makefile.llamp (modified) (2 diffs)
-
module_definitions.f90 (modified) (1 diff)
-
module_generic.f90 (modified) (3 diffs)
-
module_scientific.f90 (modified) (7 diffs)
-
trajectories_overlap.f90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/tools/Makefile.llamp
r1655 r1660 14 14 LIB_INC = 15 15 RM = rm -f 16 #DBGFLAGS = -g -Wall -Wextra -Warray-temporaries -Wconversion -fimplicit-none -fbacktrace -ffree-line-length-0 -fcheck=all -ffpe-trap=zero,overflow,underflow -finit-real=nan16 DBGFLAGS = -g -Wall -Wextra -Warray-temporaries -Wconversion -fimplicit-none -fbacktrace -ffree-line-length-0 -fcheck=all -ffpe-trap=zero,overflow,underflow -finit-real=nan 17 17 NCLIBFOLD = /usr/lib/x86_64-linux-gnu 18 NCINCFOLD = /usr 19 LIB_NETCDF = -L$(NCLIBFOLD) /lib -lnetcdff -lnetcdf -I$(NCINCFOLD)/include18 NCINCFOLD = /usr/include 19 LIB_NETCDF = -L$(NCLIBFOLD) -lnetcdff -lnetcdf -I$(NCINCFOLD) 20 20 21 21 FCFLAGS = $(FCF) $(DBGFLAGS) … … 89 89 90 90 pydistrimods.o: 91 f2py -c $(NCINCFOLD) -m module_ForDistriCorrect $(distrisrcs)91 f2py -c -I$(NCINCFOLD) -m module_ForDistriCorrect $(distrisrcs) -L$(NCLIBFOLD) 92 92 93 93 pydiagmods.o: 94 f2py -c $(NCINCFOLD) -m module_ForDiag $(diagsrcs)94 f2py -c -I$(NCINCFOLD) -m module_ForDiag $(diagsrcs) -L$(NCLIBFOLD) 95 95 96 96 pyintmods.o: 97 f2py -c $(NCINCFOLD) -m module_ForInt $(intsrcs)97 f2py -c -I$(NCINCFOLD) -m module_ForInt $(intsrcs) -L$(NCLIBFOLD) 98 98 99 99 pyscimods.o: 100 f2py -c $(NCINCFOLD) -m module_ForSci $(scisrcs)100 f2py -c -I$(NCINCFOLD) -m module_ForSci $(scisrcs) -L$(NCLIBFOLD) 101 101 102 102 trajectories_overlap.o: module_definitions.o module_basic.o module_generic.o module_scientific.o -
trunk/tools/module_definitions.f90
r1655 r1660 19 19 CHARACTER(len=50) :: fname 20 20 CHARACTER(len=256) :: msg 21 ! Some all-purpose variables 22 CHARACTER(len=1) :: Str1 23 CHARACTER(len=4) :: numSa, numSb 24 CHARACTER(len=10) :: Str10 21 25 22 26 END MODULE module_definitions -
trunk/tools/module_generic.f90
r1658 r1660 11 11 ! Index2DArrayR: Function to provide the first index of a given value inside a 2D real array 12 12 ! Index2DArrayR_K: Function to provide the first index of a given value inside a 2D real(r_k) array 13 ! Nvalues_2DArrayI: Number of different values of a 2D integer array 13 14 ! mat2DPosition: Function to provide the i, j indices of a given value inside a 2D matrix 14 15 ! RangeI: Function to provide a range of d1 values from 'iniv' to 'endv', of integer values in a vector … … 50 51 CONTAINS 51 52 53 SUBROUTINE Nvalues_2DArrayI(dx, dy, dxy, mat2DI, Nvals, vals) 54 ! Subroutine to give the number of different values of a 2D integer array 55 56 IMPLICIT NONE 57 58 INTEGER, INTENT(in) :: dx, dy, dxy 59 INTEGER, DIMENSION(dx,dy), INTENT(in) :: mat2DI 60 INTEGER, INTENT(out) :: Nvals 61 INTEGER, DIMENSION(dxy), INTENT(out) :: vals 62 63 ! Local 64 INTEGER :: i, j, ij 65 66 !!!!!!! Variables 67 ! dx, dy: size of the 2D space 68 ! mat2DI: 2D integer matrix 69 ! Nvals: number of different values 70 ! vals: vector with the different values 71 72 fname = 'Nvalues_2DArrayI' 73 74 vals = 0 75 76 Nvals = 1 77 vals(1) = mat2DI(1,1) 78 DO i=1,dx 79 DO j=1,dy 80 IF (Index1DArrayI(vals, Nvals, mat2DI(i,j)) == -1) THEN 81 Nvals = Nvals + 1 82 vals(Nvals) = mat2DI(i,j) 83 END IF 84 END DO 85 END DO 86 87 RETURN 88 89 END SUBROUTINE Nvalues_2DArrayI 90 52 91 INTEGER FUNCTION index_list_coordsI(Ncoords, coords, icoord) 53 92 ! Function to provide the index of a given coordinate within a list of integer coordinates … … 290 329 IMPLICIT NONE 291 330 292 INTEGER :: funit293 331 LOGICAL :: is_used 294 332 295 333 is_used = .true. 296 334 DO freeunit=10,100 297 INQUIRE(unit=f unit, opened=is_used)335 INQUIRE(unit=freeunit, opened=is_used) 298 336 IF (.not. is_used) EXIT 299 337 END DO -
trunk/tools/module_scientific.f90
r1655 r1660 28 28 ! poly_overlap_tracks_area: Subroutine to determine tracks of a series of consecutive 2D field with polygons using 29 29 ! maximum overlaping/coincidence filtrered by a minimal area 30 ! poly_overlap_tracks_area_ascii: Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum 31 ! overlaping/coincidence filtrered by a minimal area writting theoutput on an ASCII file (memory limitations) 30 32 ! quantilesR_K: Subroutine to provide the quantiles of a given set of values of type real 'r_k' 31 33 ! rand_sample: Subroutine to randomly sample a range of indices 34 ! read_finaltrack_ascii: Subroutine to read the final trajectories from an ASCII file 35 ! read_overlap_single_track_ascii: Subroutine to read the values for a given trajectory 36 ! read_overlap_polys_ascii: Subroutine to read from an ASCII file the associated polygons at a given time-step 37 ! read_overlap_tracks_ascii: Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step 32 38 ! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order. 33 39 ! StatsR_K: Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a 34 40 ! series of r_k numbers 35 41 ! SwapR_K*: Subroutine swaps the values of its two formal arguments. 42 ! write_finaltrack_ascii: Subroutine to read the final trajectories into an ASCII file 43 ! write_overlap_polys_ascii: Subroutine to write to an ASCII file the associated polygons at a given time-step 44 ! write_overlap_tracks_ascii: Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step 36 45 37 46 !!! *Functions/Subroutines to sort values adpated. The method used is usually referred to as "selection" method. … … 42 51 43 52 CONTAINS 53 54 SUBROUTINE read_finaltrack_ascii(funit, dt, itrack, ftrack) 55 ! Subroutine to read the final trajectories from an ASCII file 56 57 IMPLICIT NONE 58 59 INTEGER, INTENT(in) :: funit, dt, itrack 60 REAL(r_k), DIMENSION(4,dt), INTENT(out) :: ftrack 61 62 ! Local 63 INTEGER :: i, j, it 64 LOGICAL :: found 65 66 !!!!!!! Variables 67 ! funit: unit where to write the trajectory 68 ! dt: number of time-steps 69 ! itrack: trajectory to read the values 70 ! ftrack: values of the trajectory 71 72 fname = 'read_finaltrack_ascii' 73 74 ftrack = 0. 75 76 REWIND(funit) 77 78 it = 1 79 DO WHILE (.NOT.found) 80 81 READ(funit,10)ftrack(1,1), Str1, ((ftrack(i,j),Str1,i=2,4),Str1,j=1,dt) 82 IF (INT(ftrack(1,1)) == itrack) THEN 83 ftrack(1,2:dt) = ftrack(1,1) 84 found = .TRUE. 85 END IF 86 87 ! Just in case 88 IF (it >= dt) found = .TRUE. 89 90 END DO 91 92 RETURN 93 94 10 FORMAT(I10000000,1x,A1,1x,10000000(3(F20.10,A1),A1)) 95 96 END SUBROUTINE read_finaltrack_ascii 97 98 SUBROUTINE write_finaltrack_ascii(funit, dt, ftrack) 99 ! Subroutine to write the final trajectories into an ASCII file 100 101 IMPLICIT NONE 102 103 INTEGER, INTENT(in) :: funit, dt 104 REAL(r_k), DIMENSION(4,dt), INTENT(in) :: ftrack 105 106 ! Local 107 INTEGER :: i, j 108 109 !!!!!!! Variables 110 ! funit: unit where to write the trajectory 111 ! dt: number of time-steps 112 ! ftrack: values of the trajectory 113 114 fname = 'write_finaltrack_ascii' 115 WRITE(funit,10)INT(ftrack(1,1)), ';', ((ftrack(i,j), ',', i=2,4), ':', j=1,dt) 116 117 RETURN 118 119 10 FORMAT(I10,1x,A1,1x,10000000(3(F20.10,A1),A1)) 120 121 END SUBROUTINE write_finaltrack_ascii 122 123 SUBROUTINE read_overlap_single_track_ascii(funit, dt, Nxp, Nxtr, itrack, strack) 124 ! Subroutine to read the values for a given trajectory 125 126 IMPLICIT NONE 127 128 INTEGER, INTENT(in) :: funit, dt, Nxp, Nxtr, itrack 129 REAL(r_k), DIMENSION(5,Nxp,dt), INTENT(out) :: strack 130 131 ! Local 132 INTEGER :: i,j,k,l 133 INTEGER :: read_it, itt, it, Ntrcks 134 INTEGER, DIMENSION(Nxp) :: Npindep 135 LOGICAL :: looking 136 REAL(r_k), DIMENSION(5,Nxp,Nxtr) :: trcks 137 138 !!!!!!! Variables 139 ! funit: unit from where retrieve the values of the trajectory 140 ! dt: time-dimension 141 ! Nxp: maximum allowed number of polygons per time-step 142 ! Nxp: maximum allowed number of trajectories 143 ! itrack: trajectory Id to look for 144 ! strack: Values for the given trajectory 145 146 fname = 'read_overlap_single_track_ascii' 147 148 strack = 0. 149 150 REWIND(funit) 151 152 looking = .TRUE. 153 itt = 0 154 it = 1 155 DO WHILE (looking) 156 READ(funit,5,END=100)Str10, read_it 157 158 READ(funit,*)Ntrcks 159 DO i=1, Ntrcks 160 READ(funit,10)l, Str1, Npindep(i), Str1, ((trcks(k,j,i),Str1,k=1,5),Str1,j=1,Npindep(i)) 161 END DO 162 163 ! There is the desired trajectory at this time-step? 164 IF (ANY(INT(trcks(1,1,:)) == itrack)) THEN 165 itt = itt + 1 166 DO i=1, Ntrcks 167 IF (INT(trcks(1,1,i)) == itrack) THEN 168 DO j=1, Npindep(i) 169 strack(:,j,it) = trcks(:,j,i) 170 END DO 171 END IF 172 END DO 173 ELSE 174 ! It trajectory has already been initialized this is the end 175 IF (itt > 0) looking = .FALSE. 176 END IF 177 178 ! Just in case... ;) 179 IF (read_it >= dt) looking = .FALSE. 180 it = it + 1 181 182 IF (it > dt) looking = .FALSE. 183 184 END DO 185 186 100 CONTINUE 187 188 RETURN 189 190 5 FORMAT(A10,1x,I4) 191 10 FORMAT(I4,1x,A1,I4,1x,A1,1x,1000000(5(F20.10,A1),A1)) 192 193 END SUBROUTINE read_overlap_single_track_ascii 194 195 SUBROUTINE read_overlap_tracks_ascii(funit, tstep, Nxp, Ntrcks, trcks) 196 ! Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step 197 198 IMPLICIT NONE 199 200 INTEGER, INTENT(in) :: funit, tstep, Nxp 201 INTEGER, INTENT(out) :: Ntrcks 202 REAL(r_k), DIMENSION(5,Nxp,Nxp), INTENT(out) :: trcks 203 204 ! Local 205 INTEGER :: i, j, k, l, Npindep 206 INTEGER :: read_it 207 208 !!!!!!! Variables 209 ! funit: unit where to write the file 210 ! tstep: time-step to write the trajectories 211 ! Nxp: maximum number of polygons per time-step 212 ! Nrtcks: Number of trajectories of the given time-step 213 ! trcks: trajectories 214 215 fname = 'read_overlap_tracks_ascii' 216 217 Ntrcks = 0 218 trcks = 0 219 220 READ(funit,5)Str10, read_it 221 222 IF (read_it /= tstep) THEN 223 WRITE(numSa,'(I4)')read_it 224 WRITE(numSb,'(I4)')tstep 225 msg = 'File time-step;' // TRIM(numSa) // ' does not coincide with the one from program:' // & 226 TRIM(numSb) 227 END IF 228 229 READ(funit,*)Ntrcks 230 DO i=1, Ntrcks 231 READ(funit,10)l, Str1, Npindep, Str1, ((trcks(k,j,i),Str1,k=1,5),Str1,j=1,Npindep) 232 END DO 233 234 RETURN 235 236 5 FORMAT(A10,1x,I4) 237 10 FORMAT(I4,1x,A1,I4,1x,A1,1x,1000000(5(F20.10,A1),A1)) 238 239 END SUBROUTINE read_overlap_tracks_ascii 240 241 SUBROUTINE write_overlap_tracks_ascii(funit, tstep, Nxp, Ntrcks, trcks) 242 ! Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step 243 244 IMPLICIT NONE 245 246 INTEGER, INTENT(in) :: funit, tstep, Nxp, Ntrcks 247 REAL(r_k), DIMENSION(5,Nxp,Ntrcks) :: trcks 248 249 ! Local 250 INTEGER :: i, j, k, ii, Npindep, Nrealtracks 251 252 !!!!!!! Variables 253 ! funit: unit where to write the file 254 ! tstep: time-step to write the trajectories 255 ! Nxp: maximum number of polygons per time-step 256 ! Nrtcks: Number of trajectories of the given time-step 257 ! trcks: trajectories 258 259 fname = 'write_overlap_tracks_ascii' 260 261 WRITE(funit,5)'time-step:', tstep 262 263 ! Looking for the non-zero trajectories 264 Nrealtracks = 0 265 DO i=1, Ntrcks 266 Npindep = COUNT(trcks(2,:,i) /= zeroRK) 267 IF (Npindep /= 0) Nrealtracks = Nrealtracks + 1 268 END DO 269 WRITE(funit,*)Nrealtracks 270 271 ! Only writting the trajectories with values 272 ii = 1 273 DO i=1, Ntrcks 274 Npindep = COUNT(trcks(2,:,i) /= zeroRK) 275 IF (Npindep /= 0) THEN 276 WRITE(funit,10)ii,';', Npindep, ';', ((trcks(k,j,i),',',k=1,5),':',j=1,Npindep) 277 ii = ii + 1 278 END IF 279 END DO 280 281 RETURN 282 283 5 FORMAT(A10,1x,I4) 284 10 FORMAT(I4,1x,A1,I4,1x,A1,1x,1000000(5(F20.10,A1),A1)) 285 286 END SUBROUTINE write_overlap_tracks_ascii 287 288 SUBROUTINE read_overlap_polys_ascii(funit, tstep, Nxp, Nindep, SpIndep, NpIndep, pIndep) 289 ! Subroutine to read from an ASCII file the associated polygons at a given time-step 290 291 IMPLICIT NONE 292 293 INTEGER, INTENT(in) :: funit, tstep, Nxp 294 INTEGER, INTENT(out) :: Nindep 295 INTEGER, DIMENSION(Nxp), INTENT(out) :: SpIndep, NpIndep 296 INTEGER, DIMENSION(Nxp,Nxp), INTENT(out) :: pIndep 297 298 ! Local 299 INTEGER :: i, j, k 300 INTEGER :: read_it 301 302 !!!!!!! Variables 303 ! funit: unit associated to the file 304 ! tstep: time-step of the values 305 ! Nxp: allowed maximum numbe of polygons per time-step 306 ! Nindpe: Number of independent polygons at this time-step 307 ! SpIndep: Associated polygon to the independent one from the previous time-step 308 ! NpIndep: Number of associated polygons to the independent time-step 309 ! pIndep: polygons associated to a given independent polygon 310 311 fname = 'read_overlap_polys_ascii' 312 313 Nindep = 0 314 SpIndep = 0 315 NpIndep = 0 316 317 READ(funit,5)Str10, read_it 318 319 IF (read_it /= tstep) THEN 320 WRITE(numSa,'(I4)')read_it 321 WRITE(numSb,'(I4)')tstep 322 msg = 'File time-step;' // TRIM(numSa) // ' does not coincide with the one from program:' // & 323 TRIM(numSb) 324 END IF 325 326 READ(funit,*)Nindep 327 DO i=1, Nindep 328 READ(funit,10) k, Str1, SpIndep(i), Str1, NpIndep(i), Str1, (pIndep(i,j), Str1, j=1,NpIndep(i)) 329 END DO 330 331 RETURN 332 333 5 FORMAT(A10,1x,I4) 334 10 FORMAT(I4,1x,A1,1x,I4,1x,A1,1x,I4,A1,1x,100000(I4,A1)) 335 336 END SUBROUTINE read_overlap_polys_ascii 337 338 SUBROUTINE write_overlap_polys_ascii(funit, tstep, Nxp, Nindep, SpIndep, NpIndep, pIndep) 339 ! Subroutine to write into an ASCII file the associated polygons at a given time-step 340 341 IMPLICIT NONE 342 343 INTEGER, INTENT(in) :: funit, tstep, Nxp, Nindep 344 INTEGER, DIMENSION(Nindep), INTENT(in) :: SpIndep, NpIndep 345 INTEGER, DIMENSION(Nindep,Nxp), INTENT(in) :: pIndep 346 347 ! Local 348 INTEGER :: i, j 349 350 !!!!!!! Variables 351 ! funit: unit associated to the file 352 ! tstep: time-step of the values 353 ! Nxp: allowed maximum numbe of polygons per time-step 354 ! Nindpe: Number of independent polygons at this time-step 355 ! SpIndep: Associated polygon to the independent one from the previous time-step 356 ! NpIndep: Number of associated polygons to the independent time-step 357 ! pIndep: polygons associated to a given independent polygon 358 359 fname = 'write_overlap_polys_ascii' 360 361 WRITE(funit,5)'time-step:', tstep 362 WRITE(funit,*)Nindep, ' ! Number of independent polygons' 363 DO i=1, Nindep 364 WRITE(funit,10) i, ';', SpIndep(i), ';', NpIndep(i), ':', (pIndep(i,j), ',', j=1,NpIndep(i)) 365 END DO 366 367 RETURN 368 369 5 FORMAT(A10,1x,I4) 370 10 FORMAT(I4,1x,A1,1x,I4,1x,A1,1x,I4,A1,1x,100000(I4,A1)) 371 372 END SUBROUTINE write_overlap_polys_ascii 373 374 SUBROUTINE poly_overlap_tracks_area_ascii(dbg, compute, dx, dy, dt, minarea, inNallpolys, allpolys, & 375 ctrpolys, areapolys, Nmaxpoly, Nmaxtracks) 376 ! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum 377 ! overlaping/coincidence filtrered by a minimal area writting theoutput on an ASCII file (memory limitations) 378 379 IMPLICIT NONE 380 381 LOGICAL, INTENT(in) :: dbg 382 CHARACTER(LEN=*), INTENT(in) :: compute 383 INTEGER, INTENT(in) :: dx, dy, dt, Nmaxpoly, Nmaxtracks 384 INTEGER, DIMENSION(dt), INTENT(in) :: inNallpolys 385 INTEGER, DIMENSION(dx,dy,dt), INTENT(in) :: allpolys 386 REAL(r_k), INTENT(in) :: minarea 387 REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in) :: ctrpolys 388 REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in) :: areapolys 389 390 ! Local 391 INTEGER :: i, j, ip, it, iip, itt, iit 392 INTEGER :: fprevunit, ftrackunit, ftrunit, ierr, ios 393 LOGICAL :: file_exist, dooverlap, dotracks, doftracks 394 REAL(r_k), DIMENSION(Nmaxpoly) :: Aprevpolys, Acurrpolys 395 REAL(r_k), DIMENSION(2,Nmaxpoly) :: Cprevpolys, Ccurrpolys 396 INTEGER, DIMENSION(dx,dy) :: meetpolys, searchpolys 397 INTEGER, DIMENSION(Nmaxpoly) :: coincidencies 398 INTEGER, DIMENSION(Nmaxpoly) :: prevID, currID 399 REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,2) :: tracks 400 REAL(r_k), DIMENSION(4,dt) :: finaltracks 401 INTEGER, DIMENSION(:), ALLOCATABLE :: coins 402 INTEGER, DIMENSION(:,:), ALLOCATABLE :: coinsNpts 403 INTEGER :: Nmeet, Nsearch, Nindep 404 INTEGER, DIMENSION(2) :: Nindeppolys, Npolystime 405 CHARACTER(len=5) :: NcoinS 406 INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,2) :: polysIndep 407 INTEGER, DIMENSION(Nmaxpoly,2) :: NpolysIndep 408 INTEGER, DIMENSION(Nmaxpoly,2) :: SpolysIndep 409 INTEGER :: iindep, iiprev 410 INTEGER :: Nprev, NNprev, Ntprev 411 LOGICAL :: Indeppolychained 412 INTEGER :: itrack, ictrack 413 INTEGER :: ixp, iyp, ttrack 414 INTEGER, DIMENSION(2) :: Ntracks 415 INTEGER :: idtrack, maxtrack 416 REAL(r_k), DIMENSION(5,Nmaxpoly,dt) :: singletrack 417 418 !!!!!!! Variables 419 ! dx,dy,dt: space/time dimensions 420 ! compute: how to copmute 421 ! 'scratch': everything from the beginning 422 ! 'continue': skipt that parts which already have the ascii file written 423 ! inNallpolys: Vector with the original number of polygons at each time-step 424 ! allpolys: Series of 2D field with the polygons 425 ! minarea: minimal area (in same units as areapolys) to perform the tracking 426 ! ctrpolys: center of the polygons 427 ! areapolys: area of the polygons 428 ! Nmaxpoly: Maximum possible number of polygons 429 ! Nmaxtracks: maximum number of tracks 430 431 fname = 'poly_overlap_tracks_area_ascii' 432 433 IF (dbg) PRINT *,TRIM(fname) 434 435 SELECT CASE (TRIM(compute)) 436 CASE ('scratch') 437 dooverlap = .TRUE. 438 dotracks = .TRUE. 439 doftracks = .TRUE. 440 CASE ('continue') 441 INQUIRE(file='polygons_overlap.dat', exist=file_exist) 442 IF (.NOT.file_exist) THEN 443 dooverlap = .TRUE. 444 ELSE 445 IF (dbg) THEN 446 PRINT *, TRIM(warnmsg) 447 PRINT *," "//TRIM(fname) // ": File 'polygons_overlap.dat' already exists, skipping it !!" 448 END IF 449 dooverlap = .FALSE. 450 END IF 451 INQUIRE(file='trajectories_overlap.dat', exist=file_exist) 452 IF (.NOT.file_exist) THEN 453 dotracks = .TRUE. 454 ELSE 455 IF (dbg) THEN 456 PRINT *, TRIM(warnmsg) 457 PRINT *, " " // TRIM(fname) // ": File 'trajectories_overlap.dat' already exists, " // & 458 "skipping it !!" 459 END IF 460 dotracks = .FALSE. 461 END IF 462 INQUIRE(file='trajectories.dat', exist=file_exist) 463 IF (.NOT.file_exist) THEN 464 doftracks = .TRUE. 465 ELSE 466 IF (dbg) THEN 467 PRINT *, TRIM(warnmsg) 468 PRINT *," "//TRIM(fname) // ": File 'trajectories.dat' already exists, skipping it !!" 469 END IF 470 doftracks = .FALSE. 471 END IF 472 CASE DEFAULT 473 msg = "compute case: '" // TRIM(compute) // "' not ready !!" 474 CALL ErrMsg(msg, fname, -1) 475 END SELECT 476 477 IF (dooverlap) THEN 478 ! ASCII file for all the polygons and their previous associated one 479 fprevunit = freeunit() 480 OPEN(fprevunit, file='polygons_overlap.dat', status='new', form='formatted', iostat=ios) 481 msg = "Problems opening file: 'polygons_overlap.dat'" 482 IF (ios == 17) PRINT *," Careful: 'polygons_overlap.dat' already exists!!" 483 CALL ErrMsg(msg, fname, ios) 484 485 ! Number of independent polygons by time step 486 Nindeppolys = 0 487 ! Number of polygons attached to each independent polygons by time step 488 NpolysIndep = 0 489 ! ID of searching polygon attached to each independent polygons by time step 490 SpolysIndep = 0 491 ! ID of polygons attached to each independent polygons by time step 492 polysIndep = 0 493 ! ID of polygons from previous time-step 494 prevID = 0 495 ! ID of polygons from current time-step 496 currID = 0 497 498 ! First time-step all are independent polygons 499 it = 1 500 Nmeet = inNallpolys(it) 501 Nindeppolys(it) = Nmeet 502 ip = 0 503 meetpolys = allpolys(:,:,it) 504 DO i=1, Nmeet 505 IF (areapolys(i,it) >= minarea) THEN 506 ip = ip + 1 507 SpolysIndep(ip,it) = i 508 currID(ip) = i 509 Acurrpolys(ip) = areapolys(i,it) 510 Ccurrpolys(1,ip) = ctrpolys(1,i,it) 511 Ccurrpolys(2,ip) = ctrpolys(2,i,it) 512 NpolysIndep(ip,it) = 1 513 polysIndep(ip,1,it) = i 514 ELSE 515 WHERE (meetpolys == i) 516 meetpolys = 0 517 END WHERE 518 END IF 519 END DO 520 Nindeppolys(1) = ip 521 Npolystime(1) = ip 522 523 ! Starting step 524 it = 0 525 IF (dbg) THEN 526 PRINT *,' time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',0 527 PRINT *,' number of independent polygons:', Nindeppolys(it+1) 528 PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______' 529 DO i=1, Nindeppolys(it+1) 530 PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', & 531 polysIndep(i,1:NpolysIndep(i,it+1),it+1) 532 END DO 533 END IF 534 ! Writting to the ASCII file Independent polygons and their associated 535 CALL write_overlap_polys_ascii(fprevunit,it+1, Nmaxpoly, Nindeppolys(it+1), & 536 SpolysIndep(1:Nindeppolys(it+1),it+1), NpolysIndep(1:Nindeppolys(it+1),it+1), & 537 polysIndep(1:Nindeppolys(it+1),:,it+1)) 538 539 it = 1 540 ! Looking for the coincidencies at each time step 541 DO iit=1, dt-1 542 ! Number of times that a polygon has a coincidence 543 coincidencies = 0 544 545 ! Preparing for next time-step 546 searchpolys = meetpolys 547 prevID = 0 548 prevID = currID 549 Aprevpolys = Acurrpolys 550 Cprevpolys = Ccurrpolys 551 552 Nmeet = inNallpolys(iit+1) 553 meetpolys = allpolys(:,:,iit+1) 554 ip = 0 555 DO i=1, Nmeet 556 IF (areapolys(i,iit+1) >= minarea) THEN 557 ip = ip + 1 558 currID(ip) = i 559 Acurrpolys(ip) = areapolys(i,iit+1) 560 Acurrpolys(ip) = areapolys(i,iit+1) 561 Ccurrpolys(1,ip) = ctrpolys(1,i,iit+1) 562 Ccurrpolys(2,ip) = ctrpolys(2,i,iit+1) 563 ELSE 564 WHERE (meetpolys == i) 565 meetpolys = 0 566 END WHERE 567 END IF 568 END DO 569 Nindeppolys(it+1) = ip 570 Npolystime(it+1) = ip 571 572 ! Looking throughout the independent polygons 573 Nmeet = Nindeppolys(it+1) 574 !Nsearch = Nindeppolys(it) 575 ! Previous space might have more polygons that their number of independent ones 576 Nsearch = Npolystime(it) 577 578 IF (ALLOCATED(coins)) DEALLOCATE(coins) 579 ALLOCATE(coins(Nmeet), STAT=ierr) 580 msg="Problems allocating 'coins'" 581 CALL ErrMsg(msg,fname,ierr) 582 583 IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts) 584 ALLOCATE(coinsNpts(Nmeet, Nsearch), STAT=ierr) 585 msg="Problems allocating 'coinsNpts'" 586 CALL ErrMsg(msg,fname,ierr) 587 588 PRINT *,' Lluis max(searchpolys):', MAXVAL(searchpolys), ' Nsearch:', Nsearch 589 590 CALL coincidence_all_polys_area(dbg, dx, dy, Nmeet, currID, meetpolys, Ccurrpolys(:,1:Nmeet), & 591 Nsearch, prevID, searchpolys, Cprevpolys(:,1:Nsearch), Aprevpolys(1:Nsearch), coins, & 592 coinsNpts) 593 594 ! Counting the number of times that a polygon has a coincidency 595 IF (dbg) THEN 596 PRINT *,' Coincidencies for the given time-step:', iit+1,' _______' 597 DO i=1, Nmeet 598 PRINT *,currID(i), coins(i),' N search pts:', coinsNpts(i,:) 599 END DO 600 END IF 601 602 ! Looking for the same equivalencies 603 Nindep = 0 604 DO i=1, Nmeet 605 IF (coins(i) == -1) THEN 606 Nindep = Nindep + 1 607 SpolysIndep(Nindep,it+1) = -1 608 NpolysIndep(Nindep,it+1) = NpolysIndep(Nindep,it+1) + 1 609 polysIndep(Nindep,NpolysIndep(Nindep,it+1),it+1) = currID(i) 610 ELSE IF (coins(i) == -9) THEN 611 WRITE(NcoinS,'(I5)')coins(i) 612 msg="coins= "//TRIM(NcoinS)//" This is an error. One should have always only one " // & 613 "coincidence of polygon" 614 CALL ErrMsg(msg, fname, -1) 615 ELSE 616 ! Looking for coincidences with previous independent polygons 617 DO ip=1, Nsearch 618 ! Looking into the polygons associated 619 NNprev = NpolysIndep(ip,it) 620 DO j=1, NNprev 621 IF (coins(i) == polysIndep(ip,j,it)) THEN 622 ! Which index corresponds to this coincidence? 623 iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i)) 624 IF (iindep == -1) THEN 625 Nindep = Nindep + 1 626 SpolysIndep(Nindep,it+1) = coins(i) 627 END IF 628 iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i)) 629 IF (iindep < 0) THEN 630 PRINT *,' Looking for:', coins(i) 631 PRINT *,' Within:', SpolysIndep(1:Nindep,it+1) 632 PRINT *,' Might content:', polysIndep(ip,1:NNprev,it) 633 PRINT *,' From an initial list:', coins(1:Nmeet) 634 msg = 'Wrong index! There must be an index here' 635 CALL ErrMsg(msg,fname,iindep) 636 END IF 637 coincidencies(ip) = coincidencies(ip) + 1 638 NpolysIndep(iindep,it+1) = NpolysIndep(iindep,it+1) + 1 639 polysIndep(iindep,NpolysIndep(iindep,it+1),it+1) = currID(i) 640 EXIT 641 END IF 642 END DO 643 END DO 644 END IF 645 END DO 646 Nindeppolys(it+1) = Nindep 647 648 IF (dbg) THEN 649 PRINT *,' time step:',iit+1,' number to look polygons:', Nmeet,' searching polygons:',Nsearch 650 PRINT *,' number of independent polygons:', Nindeppolys(it+1) 651 PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______' 652 DO i=1, Nindeppolys(it+1) 653 PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', & 654 polysIndep(i,1:NpolysIndep(i,it+1),it+1) 655 END DO 656 END IF 657 658 ! Writting to the ASCII file Independent polygons and their associated 659 CALL write_overlap_polys_ascii(fprevunit, iit+1, Nmaxpoly, Nindeppolys(it+1), & 660 SpolysIndep(1:Nindeppolys(it+1),it+1), NpolysIndep(1:Nindeppolys(it+1),it+1), & 661 polysIndep(1:Nindeppolys(it+1),:,it+1)) 662 ! Preparing for the next time-step 663 SpolysIndep(:,it) = 0 664 NpolysIndep(:,it) = 0 665 polysIndep(:,:,it) = 0 666 Nindeppolys(it) = Nindeppolys(it+1) 667 SpolysIndep(1:Nindeppolys(it),it) = SpolysIndep(1:Nindeppolys(it+1),it+1) 668 NpolysIndep(1:Nindeppolys(it),it) = NpolysIndep(1:Nindeppolys(it+1),it+1) 669 Npolystime(it) = Npolystime(it+1) 670 671 DO ip=1, Nindeppolys(it) 672 polysIndep(ip,1,it) = polysIndep(ip,1,it+1) 673 polysIndep(ip,2,it) = polysIndep(ip,2,it+1) 674 END DO 675 SpolysIndep(:,it+1) = 0 676 NpolysIndep(:,it+1) = 0 677 polysIndep(:,:,it+1) = 0 678 679 END DO 680 CLOSE(fprevunit) 681 IF (dbg) PRINT *," Succesful writting of ASCII chain of polygons 'polygons_overlap.dat' !!" 682 END IF 683 ! ASCII file for all the polygons and their previous associated one 684 fprevunit = freeunit() 685 OPEN(fprevunit, file='polygons_overlap.dat', status='old', form='formatted', iostat=ios) 686 msg = "Problems opening file: 'polygons_overlap.dat'" 687 CALL ErrMsg(msg, fname, ios) 688 689 it = 1 690 IF (dbg) THEN 691 PRINT *, 'Coincidencies to connect _______' 692 DO iit=1, dt 693 ! Reading from the ASCII file Independent polygons and their associated 694 CALL read_overlap_polys_ascii(fprevunit, iit, Nmaxpoly, Nindeppolys(it), SpolysIndep(:,it), & 695 NpolysIndep(:,it), polysIndep(:,:,it)) 696 PRINT *,' it:', iit, ' Nindep:', Nindeppolys(it) 697 PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs' 698 DO ip=1, Nindeppolys(it) 699 PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':', & 700 polysIndep(ip,1:NpolysIndep(ip,it),it) 701 END DO 702 END DO 703 END IF 704 705 REWIND(fprevunit) 706 707 ! Trajectories 708 ! It should be done following the number of 'independent' polygons 709 ! One would concatenate that independent polygons which share IDs from one step to another 710 IF (dotracks) THEN 711 712 ! ASCII file for the trajectories 713 ftrackunit = freeunit() 714 OPEN(ftrackunit, file='trajectories_overlap.dat', status='new', form='formatted', iostat=ios) 715 msg = "Problems opening file: 'trajectories_overlap.dat'" 716 IF (ios == 17) PRINT *," Careful: 'trajectories_overlap.dat' already exists!!" 717 CALL ErrMsg(msg,fname,ios) 718 719 ! First time-step. Take all polygons 720 itrack = 0 721 tracks = 0. 722 Ntracks = 0 723 it = 1 724 iit = 1 725 CALL read_overlap_polys_ascii(fprevunit, iit, Nmaxpoly, Nindeppolys(it), SpolysIndep(:,it), & 726 NpolysIndep(:,it), polysIndep(:,:,it)) 727 728 DO ip=1, Nindeppolys(1) 729 itrack = itrack + 1 730 tracks(1,1,itrack,1) = itrack*1. 731 tracks(2,1,itrack,1) = SpolysIndep(ip,1) 732 tracks(3,1,itrack,1) = ctrpolys(1,ip,1) 733 tracks(4,1,itrack,1) = ctrpolys(2,ip,1) 734 tracks(5,1,itrack,1) = 1 735 Ntracks(1) = Ntracks(1) + 1 736 END DO 737 738 ! Writting first time-step trajectories to the intermediate file 739 CALL write_overlap_tracks_ascii(ftrackunit,iit,Nmaxpoly, Ntracks(it), tracks(:,:,1:Ntracks(it),it)) 740 741 ! Looping allover already assigned tracks 742 it = 2 743 maxtrack = Ntracks(1) 744 timesteps: DO iit=2, dt 745 CALL read_overlap_polys_ascii(fprevunit, iit, Nmaxpoly, Nindeppolys(it), SpolysIndep(:,it), & 746 NpolysIndep(:,it), polysIndep(:,:,it)) 747 IF (dbg) PRINT *,'track-timestep:', iit, 'N indep polys:', Nindeppolys(it) 748 ! Indep polygons current time-step 749 current_poly: DO i=1, Nindeppolys(it) 750 IF (dbg) PRINT *,' curent poly:', i, 'Prev poly:', SpolysIndep(i,it), ' N ass. polygons:', & 751 NpolysIndep(i,it), 'ass. poly:', polysIndep(i,1:NpolysIndep(i,it),it) 752 Indeppolychained = .FALSE. 753 754 ! Number of tracks previous time-step 755 ! Looping overall 756 it1_tracks: DO itt=1, Ntracks(it-1) 757 itrack = tracks(1,1,itt,it-1) 758 ! Number polygons ID assigned 759 Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0) 760 IF (dbg) PRINT *,itt,' track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1) 761 762 ! Looking for coincidencies 763 DO iip=1, Ntprev 764 IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN 765 Indeppolychained = .TRUE. 766 IF (dbg) PRINT *,' coincidence found by polygon:', tracks(2,iip,itt,it-1) 767 EXIT 768 END IF 769 END DO 770 IF (Indeppolychained) THEN 771 Ntracks(it) = Ntracks(it) + 1 772 ictrack = Ntracks(it) 773 ! Assigning all the IDs to the next step of the track 774 DO iip=1, NpolysIndep(i,it) 775 iiprev = polysIndep(i,iip,it) 776 tracks(1,iip,ictrack,it) = itrack 777 tracks(2,iip,ictrack,it) = iiprev 778 ixp = ctrpolys(1,iiprev,iit) 779 iyp = ctrpolys(2,iiprev,iit) 780 tracks(3,iip,ictrack,it) = ixp 781 tracks(4,iip,ictrack,it) = iyp 782 tracks(5,iip,ictrack,it) = iit 783 END DO 784 EXIT 785 END IF 786 IF (Indeppolychained) EXIT 787 END DO it1_tracks 788 789 ! Creation of a new track 790 IF (.NOT.Indeppolychained) THEN 791 Ntracks(it) = Ntracks(it) + 1 792 ictrack = Ntracks(it) 793 ! ID of new track 794 maxtrack = maxtrack + 1 795 IF (dbg) PRINT *,' New track!', maxtrack 796 797 ! Assigning all the IDs to the next step of the track 798 DO j=1, NpolysIndep(i,it) 799 iiprev = polysIndep(i,j,it) 800 tracks(1,j,ictrack,it) = maxtrack 801 tracks(2,j,ictrack,it) = iiprev 802 ixp = ctrpolys(1,iiprev,iit) 803 iyp = ctrpolys(2,iiprev,iit) 804 tracks(3,j,ictrack,it) = ixp 805 tracks(4,j,ictrack,it) = iyp 806 tracks(5,j,ictrack,it) = iit 807 END DO 808 END IF 809 810 END DO current_poly 811 812 IF (dbg) THEN 813 PRINT *,' At this time-step:', iit, ' N trajectories:', Ntracks(it) 814 DO i=1, Ntracks(it) 815 Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0) 816 PRINT *,i ,'ID tracks:', tracks(1,1,i,it), 'ID polygon:', tracks(2,1:Nprev,i,it) 817 END DO 818 END IF 819 820 CALL write_overlap_tracks_ascii(ftrackunit,iit,Nmaxpoly,Ntracks(it),tracks(:,:,1:Ntracks(it),it)) 821 ! Re-initializing for the next time-step 822 tracks(:,:,:,it-1) = 0 823 Ntracks(it-1) = Ntracks(it) 824 tracks(:,:,1:Ntracks(it-1),it-1) = tracks(:,:,1:Ntracks(it),it) 825 Ntracks(it) = 0 826 tracks(:,:,:,it) = 0 827 828 END DO timesteps 829 CLOSE(ftrackunit) 830 IF (dbg) PRINT *," Succesful writting of ASCII chain of polygons 'trajectories_overlap.dat' !!" 831 CLOSE(fprevunit) 832 END IF 833 834 ! Summarizing trajectories 835 ! When multiple polygons are available, the mean of their central positions determines the position 836 837 IF (doftracks) THEN 838 ! ASCII file for the trajectories 839 ftrackunit = freeunit() 840 OPEN(ftrackunit, file='trajectories_overlap.dat', status='old', form='formatted', iostat=ios) 841 msg = "Problems opening file: 'trajectories_overlap.dat'" 842 CALL ErrMsg(msg,fname,ios) 843 844 ! ASCII file for the final trajectories 845 ftrunit = freeunit() 846 OPEN(ftrunit, file='trajectories.dat', status='new', form='formatted', iostat=ios) 847 msg = "Problems opening file: 'trajectories.dat'" 848 IF (ios == 17) PRINT *," Careful: 'trajectories.dat' already exists!!" 849 CALL ErrMsg(msg,fname,ios) 850 851 finaltracks = 0. 852 853 DO itt=1, Nmaxtracks 854 CALL read_overlap_single_track_ascii(ftrackunit, dt, Nmaxpoly, Nmaxtracks, itt, singletrack) 855 856 ! It might reach the las trajectory 857 IF (ALL(singletrack == zeroRK)) EXIT 858 859 itrack = INT(MAXVAL(singletrack(1,1,:))) 860 IF (dbg) THEN 861 PRINT *,' Trajectory:', itt, '_______', itrack 862 DO it=1, dt 863 IF (singletrack(2,1,it) /= zeroRK) THEN 864 j = COUNT(singletrack(2,:,it) /= zeroRK) 865 PRINT *,it,':',(singletrack(3,i,it),',',singletrack(4,i,it),' ; ',i=1,j) 866 END IF 867 END DO 868 END IF 869 870 finaltracks = zeroRK 871 finaltracks(1,:) = itrack*1. 872 DO it =1, dt 873 Nprev = COUNT(INT(singletrack(2,:,it)) /= zeroRK) 874 IF (Nprev /= 0) THEN 875 finaltracks(2,it) = SUM(singletrack(3,:,it))/Nprev*1. 876 finaltracks(3,it) = SUM(singletrack(4,:,it))/Nprev*1. 877 finaltracks(4,it) = it*1. 878 END IF 879 END DO 880 ! Writting the final track into the ASCII file 881 CALL write_finaltrack_ascii(ftrunit, dt, finaltracks) 882 883 END DO 884 CLOSE(ftrackunit) 885 IF (dbg) PRINT *," Succesful writting of ASCII trajectories 'trajectories.dat' !!" 886 CLOSE(ftrunit) 887 END IF 888 889 IF (ALLOCATED(coins)) DEALLOCATE(coins) 890 IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts) 891 892 RETURN 893 894 END SUBROUTINE poly_overlap_tracks_area_ascii 44 895 45 896 SUBROUTINE poly_overlap_tracks_area(dbg, dx, dy, dt, minarea, inNallpolys, allpolys, ctrpolys, & … … 430 1281 ! Nallpoly: Number of polygons to find coincidence 431 1282 ! allpoly: space with the polygons to meet 1283 ! IDallpoly: ID of the polygon to find coincidence 432 1284 ! icpolys: center of the polygons to look for the coincidence 433 1285 ! Npoly: number of polygons on the 2D space 434 1286 ! polys: 2D field of polygons identified by their integer number (0 for no polygon) 1287 ! IDpolys: ID of the polygon to search for coincidences 435 1288 ! cpolys: center of the polygons 436 1289 ! apolys: area of the polygons … … 494 1347 INTEGER :: maxcorr 495 1348 INTEGER :: Nmaxcorr 1349 ! Lluis 1350 INTEGER :: Ndiffvals 1351 INTEGER, DIMENSION(:), ALLOCATABLE :: diffvals 496 1352 497 1353 !!!!!!! Variables … … 518 1374 PRINT *,' 2D polygons space ...' 519 1375 DO i=1,dx 520 PRINT '(1000(I 3,1x))',polys(i,:)1376 PRINT '(1000(I7,1x))',polys(i,:) 521 1377 END DO 1378 END IF 1379 1380 IF (ALLOCATED(diffvals)) DEALLOCATE(diffvals) 1381 ALLOCATE(diffvals(dx*dy)) 1382 1383 ! Checking for consistency on number of polygons and real content (except 0 value) 1384 CALL Nvalues_2DArrayI(dx, dy, dx*dy, polys, Ndiffvals, diffvals) 1385 IF (Ndiffvals -1 /= Npoly) THEN 1386 PRINT *,TRIM(emsg) 1387 PRINT *,' number of different values:', Ndiffvals-1, ' theoretical Npoly:', Npoly 1388 PRINT *,' Different values:', diffvals(1:Ndiffvals) 1389 msg = 'Number of different values and Npoly must coincide' 1390 CALL ErrMsg(msg, fname, -1) 522 1391 END IF 523 1392 … … 530 1399 IF (poly(i,j) .AND. polys(i,j) .NE. 0) THEN 531 1400 IF (.NOT.ANY(IDpoly == polys(i,j))) THEN 1401 IF ((ip) > 0) PRINT *,' Lluis ip:', ip, ' Npoly:', Npoly, ' ID:', polys(i,j), ' IDpoly:', IDpoly(1:ip), & 1402 ' coinNpts:', coinNpts(ip) 532 1403 ip = ip + 1 533 1404 IDpoly(ip) = polys(i,j) … … 964 1835 PRINT *,' 2D polygons space ...' 965 1836 DO i=1,dx 966 PRINT '(1000(I 3,1x))',polys(i,:)1837 PRINT '(1000(I7,1x))',polys(i,:) 967 1838 END DO 968 1839 END IF -
trunk/tools/trajectories_overlap.f90
r1659 r1660 17 17 18 18 ! namelist 19 INTEGER :: missing_value 19 20 REAL(r_k) :: minarea 20 21 CHARACTER(len=50) :: polyfilename 21 22 CHARACTER(len=50) :: poly2Dn, polyNn, polywctrn, polywarean, & 22 varnresx, varnresy, timen, projectn 23 varnresx, varnresy, timen, projectn, waycomp 23 24 LOGICAL :: debug 24 25 … … 34 35 REAL(r_k) :: resx, resy 35 36 REAL(r_k), DIMENSION(:,:), ALLOCATABLE :: polyas 36 REAL(r_k), DIMENSION(:,:,:), ALLOCATABLE :: polywctrs, finaltracks 37 REAL(r_k), DIMENSION(:,:,:,:), ALLOCATABLE :: polytracks 37 REAL(r_k), DIMENSION(:,:,:), ALLOCATABLE :: polywctrs 38 38 LOGICAL :: gotvar 39 39 CHARACTER(len=3) :: num3S … … 42 42 43 43 NAMELIST /files/ polyfilename, poly2Dn, polyNn, polywctrn, polywarean, varnresx, varnresy, timen 44 NAMELIST /config/ mi narea, projectn, debug44 NAMELIST /config/ missing_value, minarea, projectn, waycomp, debug 45 45 46 46 !!!!!!! … … 54 54 ! latn: name of the variable with the latitudes 55 55 ! timen: name of the variable with the times 56 ! missing_value: value to determine the missing value 56 57 ! minarea: minimal area of the polygons to determine which polygons to use 57 58 ! projectn: name of the projection of the data in the file. Available ones: … … 59 60 ! 'lon/lat': for regular longitude-latitude 60 61 ! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR) 61 ! 62 ! waycomp: how to copmute when intermediate ASCII files are used 63 ! 'scratch': everything from the beginning 64 ! 'continue': skipt that parts which already have the ascii file written 62 65 63 66 fname = 'trajectories_overlap' … … 101 104 IF (debug) PRINT *," Got variable 'polys':", UBOUND(polys) 102 105 106 ! Getting rid of the mising values 107 WHERE (polys == missing_value) 108 polys = 0 109 END WHERE 110 103 111 dimx = dims3D(1) 104 112 dimy = dims3D(2) … … 194 202 ! Preparing outcome (trajectories) 195 203 Nmaxtrack = Nmaxpoly*dimt/100 196 197 IF (ALLOCATED(polytracks)) DEALLOCATE(polytracks)198 ALLOCATE(polytracks(5,Nmaxpoly,Nmaxtrack,dimt), STAT=ierr)199 msg = "Error allocating 'polytracks'"200 CALL ErrMsg(msg, fname, ierr)201 202 IF (ALLOCATED(finaltracks)) DEALLOCATE(finaltracks)203 ALLOCATE(finaltracks(4,Nmaxtrack,dimt), STAT=ierr)204 msg = "Error allocating 'finaltracks'"205 CALL ErrMsg(msg, fname, ierr)206 204 207 205 IF (debug) THEN … … 214 212 215 213 ! Computing trajectories 216 CALL poly_overlap_tracks_area (debug, dimx, dimy, dimt, minarea, polyN, polys, polywctrs, polyas,&217 Nmaxpoly, Nmaxtrack, polytracks, finaltracks)214 CALL poly_overlap_tracks_area_ascii(debug, waycomp, dimx, dimy, dimt, minarea, polyN, polys, & 215 polywctrs, polyas, Nmaxpoly, Nmaxtrack) 218 216 219 217 ! Writting trajectory files
Note: See TracChangeset
for help on using the changeset viewer.
