Context Navigation

module_scientific.f90 @ 1750

Last change on this file since 1750 was 1747, checked in by lfita, 7 years ago

Adding:

`fill3DI_2Dvec': Subroutine to fill a 3D integer matrix from a series of indices from a 2D matrix
`fill3DR_2Dvec': Subroutine to fill a 3D float matrix from a series of indices from a 2D matrix

File size: 128.6 KB

Line
1	MODULE module_scientific
2	! Module of the scientific function/subroutines
3
4	!!!!!!! Functions & subroutines
5	! all_polygons_properties: Subroutine to determine the properties of all polygons in a 2D field:
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_all_polys_area: Subtourine to determine which is the coincident polygon when a boolean polygon is provided
10	! to a map of integer polygons filtrered by a minimal area
11	! coincidence_poly: Subtourine to determine which is the coincident polygon when a boolean polygon is provided
12	! to a map of integer polygons
13	! coincidence_poly_area: Subtourine to determine which is the coincident polygon when a boolean polygon is provided
14	! to a map of integer polygons filtrered by a minimal area
15	! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
16	! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
17	! fill3DI_2Dvec: Subroutine to fill a 3D integer matrix from a series of indices from a 2D matrix
18	! fill3DR_2Dvec: Subroutine to fill a 3D float matrix from a series of indices from a 2D matrix
19	! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon
20	! following ray casting algorithm
21	! look_clockwise_borders: Subroutine to look clock-wise for a next point within a collection of borders
22	! (limits of a region)
23	! paths_border: Subroutine to search the paths of a border field.
24	! path_properties: Subroutine to determine the properties of a path
25	! polygon_properties: Subroutine to determine the properties of a polygon (as .TRUE. matrix)
26	! polygons: Subroutine to search the polygons of a border field. FORTRAN based. 1st = 1!
27	! polygons_t: Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!
28	! poly_overlap_tracks: Subroutine to determine tracks of a series of consecutive 2D field with polygons using
29	! maximum overlaping/coincidence! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k)
30	! poly_overlap_tracks_area: Subroutine to determine tracks of a series of consecutive 2D field with polygons using
31	! maximum overlaping/coincidence filtrered by a minimal area
32	! poly_overlap_tracks_area_ascii: Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum
33	! overlaping/coincidence filtrered by a minimal area writting theoutput on an ASCII file (memory limitations)
34	! quantilesR_K: Subroutine to provide the quantiles of a given set of values of type real 'r_k'
35	! rand_sample: Subroutine to randomly sample a range of indices
36	! read_finaltrack_ascii: Subroutine to read the final trajectories from an ASCII file
37	! read_overlap_single_track_ascii: Subroutine to read the values for a given trajectory
38	! read_overlap_polys_ascii: Subroutine to read from an ASCII file the associated polygons at a given time-step
39	! read_overlap_tracks_ascii: Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step
40	! reconstruct_matrix: Subroutine to reconstruct a 2D matrix from a pair of syncronized vectors with the positions on x
41	! and y coordinates
42	! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order.
43	! StatsR_K: Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a
44	! series of r_k numbers
45	! SwapR_K*: Subroutine swaps the values of its two formal arguments.
46	! write_finaltrack_ascii: Subroutine to read the final trajectories into an ASCII file
47	! write_overlap_polys_ascii: Subroutine to write to an ASCII file the associated polygons at a given time-step
48	! write_overlap_tracks_ascii: Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step
49
50	!!! *Functions/Subroutines to sort values adpated. The method used is usually referred to as "selection" method.
51	! from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
52
53	USE module_definitions
54	USE module_generic
55
56	CONTAINS
57
58	SUBROUTINE fill3DI_2Dvec(matind, inmat, id1, id2, od1, od2, od3, outmat)
59	! Subroutine to fill a 3D integer matrix from a series of indices from a 2D matrix
60
61	IMPLICIT NONE
62
63	INTEGER, INTENT(in) :: id1, id2, od1, od2, od3
64	INTEGER, DIMENSION(id1,id2), INTENT(in) :: inmat
65	INTEGER, DIMENSION(od1,od2), INTENT(in) :: matind
66	INTEGER, DIMENSION(od1,od2,od3), INTENT(out) :: outmat
67
68	! Local
69	INTEGER :: i, j
70
71	!!!!!!! Variables
72	! matind: equivalence on the 2D space of the location in inmat
73	! inmat: values of the input matrix (1Dvec, time)
74	! id1/2: dimensions of the input matrix
75	! od1/3: dimensions of the output matrix
76	! outmat: output matrix
77	! NOTE: id2 == od3
78
79	fname = 'fill3DR_2Dvec'
80
81	DO i=1, od1
82	DO j=1, od2
83	IF (matind(i,j) /= -1) THEN
84	outmat(i,j,:) = inmat(matind(i,j),:)
85	ELSE
86	outmat(i,j,:) = fillvalI
87	END IF
88	END DO
89	END DO
90
91	END SUBROUTINE fill3DI_2Dvec
92
93	SUBROUTINE fill3DR_2Dvec(matind, inmat, id1, id2, od1, od2, od3, outmat)
94	! Subroutine to fill a 3D float matrix from a series of indices from a 2D matrix
95
96	IMPLICIT NONE
97
98	INTEGER, INTENT(in) :: id1, id2, od1, od2, od3
99	REAL(r_k), DIMENSION(id1,id2), INTENT(in) :: inmat
100	INTEGER, DIMENSION(od1,od2), INTENT(in) :: matind
101	REAL(r_k), DIMENSION(od1,od2,od3), INTENT(out) :: outmat
102
103	! Local
104	INTEGER :: i, j
105
106	!!!!!!! Variables
107	! matind: equivalence on the 2D space of the location in inmat
108	! inmat: values of the input matrix (1Dvec, time)
109	! id1/2: dimensions of the input matrix
110	! od1/3: dimensions of the output matrix
111	! outmat: output matrix
112	! NOTE: id2 == od3
113
114	fname = 'fill3DR_2Dvec'
115
116	DO i=1, od1
117	DO j=1, od2
118	IF (matind(i,j) /= -1) THEN
119	outmat(i,j,:) = inmat(matind(i,j),:)
120	ELSE
121	outmat(i,j,:) = fillval64
122	END IF
123	END DO
124	END DO
125
126	END SUBROUTINE fill3DR_2Dvec
127
128	SUBROUTINE reconstruct_matrix(vectorXpos, vectorYpos, dvec, Xmin, Xmax, Ymin, Ymax, dmatx, dmaty, &
129	matProj, maxdiff, matind, matX, matY)
130	! Subroutine to reconstruct a 2D matrix from a pair of syncronized vectors with the positions on x
131	! and y coordinates
132
133	IMPLICIT NONE
134
135	INTEGER, INTENT(in) :: dvec, dmatx, dmaty
136	REAL(r_k), DIMENSION(dvec), INTENT(in) :: vectorXpos, vectorYpos
137	REAL(r_k), INTENT(in) :: Xmin, Xmax, Ymin, Ymax, maxdiff
138	CHARACTER(len=50), INTENT(in) :: matPRoj
139	INTEGER, DIMENSION(dmatx, dmaty), INTENT(out) :: matind
140	REAL(r_k), DIMENSION(dmatx, dmaty), INTENT(out) :: matX, matY
141
142	! Local
143	INTEGER :: i,j, idiff, jdiff
144	REAL(r_k) :: ddx, ddy, Xpos, Ypos, mindiff
145	REAL(r_k), DIMENSION(dvec) :: diff
146	INTEGER, DIMENSION(1) :: mindiffloc
147	CHARACTER(LEN=50) :: RSa, RSb
148
149	!!!!!!! Variables
150	! vectorXpos, vectorYpos: syncronized vectors with the x,y coordinates from which the matrix will be reconstructed
151	! dvec: quantitiy of coordinates to use
152	! Xmin,Xmax,Ymin,Ymax: Location range of the matrix to be reconstructed
153	! maxdiff: Authorized maximum distance between matrix location and vector location
154	! dmatx, dmaty: Size in grid points of the matrix to be reconstructed
155	! matProj: Assumed projection of the values of the matrix
156	! 'latlon': regular lat/lon projection
157	! matind: matrix with the correspondant indiced from the vector of positions
158	! matX, matY: matrix with the coordinates of the elements of the matrix
159
160	fname = 'reconstruct_matrix'
161
162	matind = zeroRK
163
164	Projection: SELECT CASE(TRIM(matProj))
165
166	CASE('latlon')
167	! Resolution along matrix axes
168	ddx = (Xmax - Xmin)/(dmatx-1)
169	ddy = (Ymax - Ymin)/(dmaty-1)
170
171	! Filling matrix positions
172	DO i=1,dmatx
173	Xpos = Xmin + ddx*(i-1)
174	DO j=1,dmaty
175	Ypos = Ymin + ddy*(j-1)
176	matX(i,j) = Xpos
177	matY(i,j) = Ypos
178
179	diff = SQRT((Xpos - vectorXpos)2 + (Ypos - vectorYpos)2)
180	mindiffloc = MINLOC(diff)
181	mindiff = MINVAL(diff)
182	IF (mindiff > maxdiff) THEN
183	!Assuming no values
184	matind(i,j) = -1
185	!WRITE(RSa, '(f20.10)')mindiff
186	!WRITE(RSb, '(f20.10)')maxdiff
187	!msg = 'Difference: ' // TRIM(RSa) // ' overpasses the maximum authorized distance: ' // &
188	! TRIM(RSb)
189	!CALL ErrMsg(msg, fname, -1)
190	ELSE
191	matind(i,j) = mindiffloc(1)
192	END IF
193
194	END DO
195	END DO
196
197	CASE DEFAULT
198	msg = "Projection of matrix '" // TRIM(matProj) // "' not ready " // CHAR(10) // &
199	" Available ones: 'latlon'"
200	CALL ErrMsg(msg, fname, -1)
201	END SELECT Projection
202
203	END SUBROUTINE reconstruct_matrix
204
205	SUBROUTINE read_finaltrack_ascii(funit, dt, itrack, ftrack)
206	! Subroutine to read the final trajectories from an ASCII file
207
208	IMPLICIT NONE
209
210	INTEGER, INTENT(in) :: funit, dt, itrack
211	REAL(r_k), DIMENSION(5,dt), INTENT(out) :: ftrack
212
213	! Local
214	INTEGER :: i, j, it
215	LOGICAL :: found
216
217	!!!!!!! Variables
218	! funit: unit where to write the trajectory
219	! dt: number of time-steps
220	! itrack: trajectory to read the values
221	! ftrack: values of the trajectory
222
223	fname = 'read_finaltrack_ascii'
224
225	ftrack = 0.
226
227	REWIND(funit)
228
229	it = 1
230	DO WHILE (.NOT.found)
231
232	READ(funit,10)ftrack(1,1), Str1, ((ftrack(i,j),Str1,i=2,5),Str1,j=1,dt)
233	IF (INT(ftrack(1,1)) == itrack) THEN
234	ftrack(1,2:dt) = ftrack(1,1)
235	found = .TRUE.
236	END IF
237
238	! Just in case
239	IF (it >= dt) found = .TRUE.
240
241	END DO
242
243	RETURN
244
245	10 FORMAT(I10000000,1x,A1,1x,10000000(4(F20.10,A1),A1))
246
247	END SUBROUTINE read_finaltrack_ascii
248
249	SUBROUTINE write_finaltrack_ascii(funit, dt, ftrack)
250	! Subroutine to write the final trajectories into an ASCII file
251
252	IMPLICIT NONE
253
254	INTEGER, INTENT(in) :: funit, dt
255	REAL(r_k), DIMENSION(5,dt), INTENT(in) :: ftrack
256
257	! Local
258	INTEGER :: i, j
259
260	!!!!!!! Variables
261	! funit: unit where to write the trajectory
262	! dt: number of time-steps
263	! ftrack: values of the trajectory
264
265	fname = 'write_finaltrack_ascii'
266	WRITE(funit,10)INT(ftrack(1,1)), ';', ((ftrack(i,j), ',', i=2,5), ':', j=1,dt)
267
268	RETURN
269
270	10 FORMAT(I10,1x,A1,1x,10000000(4(F20.10,A1),A1))
271
272	END SUBROUTINE write_finaltrack_ascii
273
274	SUBROUTINE read_overlap_single_track_ascii(funit, dt, Nxp, Nxtr, itrack, strack)
275	! Subroutine to read the values for a given trajectory
276
277	IMPLICIT NONE
278
279	INTEGER, INTENT(in) :: funit, dt, Nxp, Nxtr, itrack
280	REAL(r_k), DIMENSION(5,Nxp,dt), INTENT(out) :: strack
281
282	! Local
283	INTEGER :: i,j,k,l
284	INTEGER :: read_it, itt, it, Ntrcks
285	INTEGER, DIMENSION(Nxp) :: Npindep
286	LOGICAL :: looking
287	REAL(r_k), DIMENSION(5,Nxp,Nxtr) :: trcks
288
289	!!!!!!! Variables
290	! funit: unit from where retrieve the values of the trajectory
291	! dt: time-dimension
292	! Nxp: maximum allowed number of polygons per time-step
293	! Nxp: maximum allowed number of trajectories
294	! itrack: trajectory Id to look for
295	! strack: Values for the given trajectory
296
297	fname = 'read_overlap_single_track_ascii'
298
299	strack = 0.
300
301	REWIND(funit)
302
303	looking = .TRUE.
304	itt = 0
305	it = 1
306	DO WHILE (looking)
307	READ(funit,5,END=100)Str10, read_it
308
309	READ(funit,*)Ntrcks
310	DO i=1, Ntrcks
311	READ(funit,10)l, Str1, Npindep(i), Str1, ((trcks(k,j,i),Str1,k=1,5),Str1,j=1,Npindep(i))
312	END DO
313
314	! There is the desired trajectory at this time-step?
315	IF (ANY(INT(trcks(1,1,:)) == itrack)) THEN
316	itt = itt + 1
317	DO i=1, Ntrcks
318	IF (INT(trcks(1,1,i)) == itrack) THEN
319	DO j=1, Npindep(i)
320	strack(:,j,it) = trcks(:,j,i)
321	END DO
322	END IF
323	END DO
324	ELSE
325	! It trajectory has already been initialized this is the end
326	IF (itt > 0) looking = .FALSE.
327	END IF
328
329	! Just in case... ;)
330	IF (read_it >= dt) looking = .FALSE.
331	it = it + 1
332
333	IF (it > dt) looking = .FALSE.
334
335	END DO
336
337	100 CONTINUE
338
339	RETURN
340
341	5 FORMAT(A10,1x,I4)
342	10 FORMAT(I4,1x,A1,I4,1x,A1,1x,1000000(5(F20.10,A1),A1))
343
344	END SUBROUTINE read_overlap_single_track_ascii
345
346	SUBROUTINE read_overlap_tracks_ascii(funit, tstep, Nxp, Ntrcks, trcks)
347	! Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step
348
349	IMPLICIT NONE
350
351	INTEGER, INTENT(in) :: funit, tstep, Nxp
352	INTEGER, INTENT(out) :: Ntrcks
353	REAL(r_k), DIMENSION(5,Nxp,Nxp), INTENT(out) :: trcks
354
355	! Local
356	INTEGER :: i, j, k, l, Npindep
357	INTEGER :: read_it
358
359	!!!!!!! Variables
360	! funit: unit where to write the file
361	! tstep: time-step to write the trajectories
362	! Nxp: maximum number of polygons per time-step
363	! Nrtcks: Number of trajectories of the given time-step
364	! trcks: trajectories
365
366	fname = 'read_overlap_tracks_ascii'
367
368	Ntrcks = 0
369	trcks = 0
370
371	READ(funit,5)Str10, read_it
372
373	IF (read_it /= tstep) THEN
374	WRITE(numSa,'(I4)')read_it
375	WRITE(numSb,'(I4)')tstep
376	msg = 'File time-step;' // TRIM(numSa) // ' does not coincide with the one from program:' // &
377	TRIM(numSb)
378	END IF
379
380	READ(funit,*)Ntrcks
381	DO i=1, Ntrcks
382	READ(funit,10)l, Str1, Npindep, Str1, ((trcks(k,j,i),Str1,k=1,5),Str1,j=1,Npindep)
383	END DO
384
385	RETURN
386
387	5 FORMAT(A10,1x,I4)
388	10 FORMAT(I4,1x,A1,I4,1x,A1,1x,1000000(5(F20.10,A1),A1))
389
390	END SUBROUTINE read_overlap_tracks_ascii
391
392	SUBROUTINE write_overlap_tracks_ascii(funit, tstep, Nxp, Ntrcks, trcks)
393	! Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step
394
395	IMPLICIT NONE
396
397	INTEGER, INTENT(in) :: funit, tstep, Nxp, Ntrcks
398	REAL(r_k), DIMENSION(5,Nxp,Ntrcks) :: trcks
399
400	! Local
401	INTEGER :: i, j, k, ii, Npindep, Nrealtracks
402
403	!!!!!!! Variables
404	! funit: unit where to write the file
405	! tstep: time-step to write the trajectories
406	! Nxp: maximum number of polygons per time-step
407	! Nrtcks: Number of trajectories of the given time-step
408	! trcks: trajectories
409
410	fname = 'write_overlap_tracks_ascii'
411
412	WRITE(funit,5)'time-step:', tstep
413
414	! Looking for the non-zero trajectories
415	Nrealtracks = 0
416	DO i=1, Ntrcks
417	Npindep = COUNT(trcks(2,:,i) /= zeroRK)
418	IF (Npindep /= 0) Nrealtracks = Nrealtracks + 1
419	END DO
420	WRITE(funit,*)Nrealtracks
421
422	! Only writting the trajectories with values
423	ii = 1
424	DO i=1, Ntrcks
425	Npindep = COUNT(trcks(2,:,i) /= zeroRK)
426	IF (Npindep /= 0) THEN
427	WRITE(funit,10)ii,';', Npindep, ';', ((trcks(k,j,i),',',k=1,5),':',j=1,Npindep)
428	ii = ii + 1
429	END IF
430	END DO
431
432	RETURN
433
434	5 FORMAT(A10,1x,I4)
435	10 FORMAT(I4,1x,A1,I4,1x,A1,1x,1000000(5(F20.10,A1),A1))
436
437	END SUBROUTINE write_overlap_tracks_ascii
438
439	SUBROUTINE read_overlap_polys_ascii(funit, tstep, Nxp, Nindep, SpIndep, NpIndep, pIndep)
440	! Subroutine to read from an ASCII file the associated polygons at a given time-step
441
442	IMPLICIT NONE
443
444	INTEGER, INTENT(in) :: funit, tstep, Nxp
445	INTEGER, INTENT(out) :: Nindep
446	INTEGER, DIMENSION(Nxp), INTENT(out) :: SpIndep, NpIndep
447	INTEGER, DIMENSION(Nxp,Nxp), INTENT(out) :: pIndep
448
449	! Local
450	INTEGER :: i, j, k
451	INTEGER :: read_it
452
453	!!!!!!! Variables
454	! funit: unit associated to the file
455	! tstep: time-step of the values
456	! Nxp: allowed maximum numbe of polygons per time-step
457	! Nindpe: Number of independent polygons at this time-step
458	! SpIndep: Associated polygon to the independent one from the previous time-step
459	! NpIndep: Number of associated polygons to the independent time-step
460	! pIndep: polygons associated to a given independent polygon
461
462	fname = 'read_overlap_polys_ascii'
463
464	Nindep = 0
465	SpIndep = 0
466	NpIndep = 0
467
468	READ(funit,5)Str10, read_it
469
470	IF (read_it /= tstep) THEN
471	WRITE(numSa,'(I4)')read_it
472	WRITE(numSb,'(I4)')tstep
473	msg = 'File time-step;' // TRIM(numSa) // ' does not coincide with the one from program:' // &
474	TRIM(numSb)
475	END IF
476
477	READ(funit,*)Nindep
478	DO i=1, Nindep
479	READ(funit,10) k, Str1, SpIndep(i), Str1, NpIndep(i), Str1, (pIndep(i,j), Str1, j=1,NpIndep(i))
480	END DO
481
482	RETURN
483
484	5 FORMAT(A10,1x,I4)
485	10 FORMAT(I4,1x,A1,1x,I4,1x,A1,1x,I4,A1,1x,100000(I4,A1))
486
487	END SUBROUTINE read_overlap_polys_ascii
488
489	SUBROUTINE write_overlap_polys_ascii(funit, tstep, Nxp, Nindep, SpIndep, NpIndep, pIndep)
490	! Subroutine to write into an ASCII file the associated polygons at a given time-step
491
492	IMPLICIT NONE
493
494	INTEGER, INTENT(in) :: funit, tstep, Nxp, Nindep
495	INTEGER, DIMENSION(Nindep), INTENT(in) :: SpIndep, NpIndep
496	INTEGER, DIMENSION(Nindep,Nxp), INTENT(in) :: pIndep
497
498	! Local
499	INTEGER :: i, j
500
501	!!!!!!! Variables
502	! funit: unit associated to the file
503	! tstep: time-step of the values
504	! Nxp: allowed maximum numbe of polygons per time-step
505	! Nindpe: Number of independent polygons at this time-step
506	! SpIndep: Associated polygon to the independent one from the previous time-step
507	! NpIndep: Number of associated polygons to the independent time-step
508	! pIndep: polygons associated to a given independent polygon
509
510	fname = 'write_overlap_polys_ascii'
511
512	WRITE(funit,5)'time-step:', tstep
513	WRITE(funit,*)Nindep, ' ! Number of independent polygons'
514	DO i=1, Nindep
515	WRITE(funit,10) i, ';', SpIndep(i), ';', NpIndep(i), ':', (pIndep(i,j), ',', j=1,NpIndep(i))
516	END DO
517
518	RETURN
519
520	5 FORMAT(A10,1x,I4)
521	10 FORMAT(I4,1x,A1,1x,I4,1x,A1,1x,I4,A1,1x,100000(I4,A1))
522
523	END SUBROUTINE write_overlap_polys_ascii
524
525	SUBROUTINE poly_overlap_tracks_area_ascii(dbg, compute, dx, dy, dt, minarea, inNallpolys, allpolys, &
526	ctrpolys, areapolys, Nmaxpoly, Nmaxtracks, methodmulti)
527	! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum
528	! overlaping/coincidence filtrered by a minimal area writting theoutput on an ASCII file (memory limitations)
529
530	IMPLICIT NONE
531
532	LOGICAL, INTENT(in) :: dbg
533	CHARACTER(LEN=*), INTENT(in) :: compute, methodmulti
534	INTEGER, INTENT(in) :: dx, dy, dt, Nmaxpoly, Nmaxtracks
535	INTEGER, DIMENSION(dt), INTENT(in) :: inNallpolys
536	INTEGER, DIMENSION(dx,dy,dt), INTENT(in) :: allpolys
537	REAL(r_k), INTENT(in) :: minarea
538	REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in) :: ctrpolys
539	REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in) :: areapolys
540
541	! Local
542	INTEGER :: i, j, ip, it, iip, itt, iit
543	INTEGER :: fprevunit, ftrackunit, ftrunit, ierr, ios
544	LOGICAL :: file_exist, dooverlap, dotracks, doftracks
545	REAL(r_k), DIMENSION(Nmaxpoly) :: Aprevpolys, Acurrpolys
546	REAL(r_k), DIMENSION(2,Nmaxpoly) :: Cprevpolys, Ccurrpolys
547	INTEGER, DIMENSION(dx,dy) :: meetpolys, searchpolys
548	INTEGER, DIMENSION(Nmaxpoly) :: coincidencies
549	INTEGER, DIMENSION(Nmaxpoly) :: prevID, currID
550	REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,2) :: tracks
551	REAL(r_k), DIMENSION(5,dt) :: finaltracks
552	INTEGER, DIMENSION(:), ALLOCATABLE :: coins
553	INTEGER, DIMENSION(:,:), ALLOCATABLE :: coinsNpts
554	INTEGER :: Nmeet, Nsearch, Nindep
555	INTEGER, DIMENSION(2) :: Nindeppolys, Npolystime
556	CHARACTER(len=5) :: NcoinS
557	INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,2) :: polysIndep
558	INTEGER, DIMENSION(Nmaxpoly,2) :: NpolysIndep
559	INTEGER, DIMENSION(Nmaxpoly,2) :: SpolysIndep
560	INTEGER :: iindep, iiprev
561	INTEGER :: Nprev, NNprev, Ntprev
562	LOGICAL :: Indeppolychained
563	INTEGER :: itrack, ictrack
564	INTEGER :: ixp, iyp, ttrack
565	INTEGER, DIMENSION(2) :: Ntracks
566	INTEGER :: idtrack, maxtrack
567	REAL(r_k), DIMENSION(5,Nmaxpoly,dt) :: singletrack
568	REAL(r_k) :: totArea, dist, mindist, maxarea, areai
569
570	!!!!!!! Variables
571	! dx,dy,dt: space/time dimensions
572	! compute: how to copmute
573	! 'scratch': everything from the beginning
574	! 'continue': skipt that parts which already have the ascii file written
575	! inNallpolys: Vector with the original number of polygons at each time-step
576	! allpolys: Series of 2D field with the polygons
577	! minarea: minimal area (in same units as areapolys) to perform the tracking
578	! ctrpolys: center of the polygons
579	! areapolys: area of the polygons
580	! Nmaxpoly: Maximum possible number of polygons
581	! Nmaxtracks: maximum number of tracks
582	! methodmulti: methodology to follow when multiple polygons are given for the same track
583	! 'mean': get coordinates from the areal-weighted mean of the centers of the given polygons and their areas
584	! 'largest': get the coorindates of the largest polygon
585	! 'closest': get the coordinates of the closest polygon
586
587	fname = 'poly_overlap_tracks_area_ascii'
588
589	IF (dbg) PRINT *,TRIM(fname)
590
591	SELECT CASE (TRIM(compute))
592	CASE ('scratch')
593	dooverlap = .TRUE.
594	dotracks = .TRUE.
595	doftracks = .TRUE.
596	CASE ('continue')
597	INQUIRE(file='polygons_overlap.dat', exist=file_exist)
598	IF (.NOT.file_exist) THEN
599	dooverlap = .TRUE.
600	ELSE
601	IF (dbg) THEN
602	PRINT *, TRIM(warnmsg)
603	PRINT *," "//TRIM(fname) // ": File 'polygons_overlap.dat' already exists, skipping it !!"
604	END IF
605	dooverlap = .FALSE.
606	END IF
607	INQUIRE(file='trajectories_overlap.dat', exist=file_exist)
608	IF (.NOT.file_exist) THEN
609	dotracks = .TRUE.
610	ELSE
611	IF (dbg) THEN
612	PRINT *, TRIM(warnmsg)
613	PRINT *, " " // TRIM(fname) // ": File 'trajectories_overlap.dat' already exists, " // &
614	"skipping it !!"
615	END IF
616	dotracks = .FALSE.
617	END IF
618	INQUIRE(file='trajectories.dat', exist=file_exist)
619	IF (.NOT.file_exist) THEN
620	doftracks = .TRUE.
621	ELSE
622	IF (dbg) THEN
623	PRINT *, TRIM(warnmsg)
624	PRINT *," "//TRIM(fname) // ": File 'trajectories.dat' already exists, skipping it !!"
625	END IF
626	doftracks = .FALSE.
627	END IF
628	CASE DEFAULT
629	msg = "compute case: '" // TRIM(compute) // "' not ready !!"
630	CALL ErrMsg(msg, fname, -1)
631	END SELECT
632
633	! Checking multi-polygon methodology
634	IF ( (TRIM(methodmulti) /= 'mean') .AND. (TRIM(methodmulti) /= 'largest') .AND. &
635	(TRIM(methodmulti) /= 'closest')) THEN
636	msg= "methodology for multiple-polygons: '"//TRIM(methodmulti)//"' not ready" // NEW_LINE('a')//&
637	" available ones: 'mean', 'largest', 'closest'"
638	CALL ErrMsg(msg, fname, -1)
639	END IF
640
641	IF (dooverlap) THEN
642	! ASCII file for all the polygons and their previous associated one
643	fprevunit = freeunit()
644	OPEN(fprevunit, file='polygons_overlap.dat', status='new', form='formatted', iostat=ios)
645	msg = "Problems opening file: 'polygons_overlap.dat'"
646	IF (ios == 17) PRINT *," Careful: 'polygons_overlap.dat' already exists!!"
647	CALL ErrMsg(msg, fname, ios)
648
649	! Number of independent polygons by time step
650	Nindeppolys = 0
651	! Number of polygons attached to each independent polygons by time step
652	NpolysIndep = 0
653	! ID of searching polygon attached to each independent polygons by time step
654	SpolysIndep = 0
655	! ID of polygons attached to each independent polygons by time step
656	polysIndep = 0
657	! ID of polygons from previous time-step
658	prevID = 0
659	! ID of polygons from current time-step
660	currID = 0
661
662	! First time-step all are independent polygons
663	it = 1
664	Nmeet = inNallpolys(it)
665	Nindeppolys(it) = Nmeet
666	ip = 0
667	meetpolys = allpolys(:,:,it)
668	DO i=1, Nmeet
669	IF (areapolys(i,it) >= minarea) THEN
670	ip = ip + 1
671	SpolysIndep(ip,it) = i
672	currID(ip) = i
673	Acurrpolys(ip) = areapolys(i,it)
674	Ccurrpolys(1,ip) = ctrpolys(1,i,it)
675	Ccurrpolys(2,ip) = ctrpolys(2,i,it)
676	NpolysIndep(ip,it) = 1
677	polysIndep(ip,1,it) = i
678	ELSE
679	WHERE (meetpolys == i)
680	meetpolys = 0
681	END WHERE
682	END IF
683	END DO
684	Nindeppolys(1) = ip
685	Npolystime(1) = ip
686
687	! Starting step
688	it = 0
689	IF (dbg) THEN
690	PRINT *,' time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',0
691	PRINT *,' number of independent polygons:', Nindeppolys(it+1)
692	PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______'
693	DO i=1, Nindeppolys(it+1)
694	PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', &
695	polysIndep(i,1:NpolysIndep(i,it+1),it+1)
696	END DO
697	END IF
698	! Writting to the ASCII file Independent polygons and their associated
699	CALL write_overlap_polys_ascii(fprevunit,it+1, Nmaxpoly, Nindeppolys(it+1), &
700	SpolysIndep(1:Nindeppolys(it+1),it+1), NpolysIndep(1:Nindeppolys(it+1),it+1), &
701	polysIndep(1:Nindeppolys(it+1),:,it+1))
702
703	it = 1
704	! Looking for the coincidencies at each time step
705	DO iit=1, dt-1
706	! Number of times that a polygon has a coincidence
707	coincidencies = 0
708
709	! Preparing for next time-step
710	searchpolys = meetpolys
711	prevID = 0
712	prevID = currID
713	Aprevpolys = Acurrpolys
714	Cprevpolys = Ccurrpolys
715
716	Nmeet = inNallpolys(iit+1)
717	meetpolys = allpolys(:,:,iit+1)
718	ip = 0
719	DO i=1, Nmeet
720	IF (areapolys(i,iit+1) >= minarea) THEN
721	ip = ip + 1
722	currID(ip) = i
723	Acurrpolys(ip) = areapolys(i,iit+1)
724	Acurrpolys(ip) = areapolys(i,iit+1)
725	Ccurrpolys(1,ip) = ctrpolys(1,i,iit+1)
726	Ccurrpolys(2,ip) = ctrpolys(2,i,iit+1)
727	ELSE
728	WHERE (meetpolys == i)
729	meetpolys = 0
730	END WHERE
731	END IF
732	END DO
733	Nindeppolys(it+1) = ip
734	Npolystime(it+1) = ip
735
736	! Looking throughout the independent polygons
737	Nmeet = Nindeppolys(it+1)
738	!Nsearch = Nindeppolys(it)
739	! Previous space might have more polygons that their number of independent ones
740	Nsearch = Npolystime(it)
741
742	IF (ALLOCATED(coins)) DEALLOCATE(coins)
743	ALLOCATE(coins(Nmeet), STAT=ierr)
744	msg="Problems allocating 'coins'"
745	CALL ErrMsg(msg,fname,ierr)
746
747	IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts)
748	ALLOCATE(coinsNpts(Nmeet, Nsearch), STAT=ierr)
749	msg="Problems allocating 'coinsNpts'"
750	CALL ErrMsg(msg,fname,ierr)
751
752	CALL coincidence_all_polys_area(dbg, dx, dy, Nmeet, currID, meetpolys, Ccurrpolys(:,1:Nmeet), &
753	Nsearch, prevID, searchpolys, Cprevpolys(:,1:Nsearch), Aprevpolys(1:Nsearch), coins, &
754	coinsNpts)
755
756	! Counting the number of times that a polygon has a coincidency
757	IF (dbg) THEN
758	PRINT *,' Coincidencies for the given time-step:', iit+1,' _______'
759	DO i=1, Nmeet
760	PRINT *,currID(i), coins(i),' N search pts:', coinsNpts(i,:)
761	END DO
762	END IF
763
764	! Looking for the same equivalencies
765	Nindep = 0
766	DO i=1, Nmeet
767	IF (coins(i) == -1) THEN
768	Nindep = Nindep + 1
769	SpolysIndep(Nindep,it+1) = -1
770	NpolysIndep(Nindep,it+1) = NpolysIndep(Nindep,it+1) + 1
771	polysIndep(Nindep,NpolysIndep(Nindep,it+1),it+1) = currID(i)
772	ELSE IF (coins(i) == -9) THEN
773	WRITE(NcoinS,'(I5)')coins(i)
774	msg="coins= "//TRIM(NcoinS)//" This is an error. One should have always only one " // &
775	"coincidence of polygon"
776	CALL ErrMsg(msg, fname, -1)
777	ELSE
778	! Looking for coincidences with previous independent polygons
779	DO ip=1, Nsearch
780	! Looking into the polygons associated
781	NNprev = NpolysIndep(ip,it)
782	DO j=1, NNprev
783	IF (coins(i) == polysIndep(ip,j,it)) THEN
784	! Which index corresponds to this coincidence?
785	iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
786	IF (iindep == -1) THEN
787	Nindep = Nindep + 1
788	SpolysIndep(Nindep,it+1) = coins(i)
789	END IF
790	iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
791	IF (iindep < 0) THEN
792	PRINT *,' Looking for:', coins(i)
793	PRINT *,' Within:', SpolysIndep(1:Nindep,it+1)
794	PRINT *,' Might content:', polysIndep(ip,1:NNprev,it)
795	PRINT *,' From an initial list:', coins(1:Nmeet)
796	msg = 'Wrong index! There must be an index here'
797	CALL ErrMsg(msg,fname,iindep)
798	END IF
799	coincidencies(ip) = coincidencies(ip) + 1
800	NpolysIndep(iindep,it+1) = NpolysIndep(iindep,it+1) + 1
801	polysIndep(iindep,NpolysIndep(iindep,it+1),it+1) = currID(i)
802	EXIT
803	END IF
804	END DO
805	END DO
806	END IF
807	END DO
808	Nindeppolys(it+1) = Nindep
809
810	IF (dbg) THEN
811	PRINT *,' time step:',iit+1,' number to look polygons:', Nmeet,' searching polygons:',Nsearch
812	PRINT *,' number of independent polygons:', Nindeppolys(it+1)
813	PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______'
814	DO i=1, Nindeppolys(it+1)
815	PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', &
816	polysIndep(i,1:NpolysIndep(i,it+1),it+1)
817	END DO
818	END IF
819
820	! Writting to the ASCII file Independent polygons and their associated
821	CALL write_overlap_polys_ascii(fprevunit, iit+1, Nmaxpoly, Nindeppolys(it+1), &
822	SpolysIndep(1:Nindeppolys(it+1),it+1), NpolysIndep(1:Nindeppolys(it+1),it+1), &
823	polysIndep(1:Nindeppolys(it+1),:,it+1))
824	! Preparing for the next time-step
825	SpolysIndep(:,it) = 0
826	NpolysIndep(:,it) = 0
827	polysIndep(:,:,it) = 0
828	Nindeppolys(it) = Nindeppolys(it+1)
829	SpolysIndep(1:Nindeppolys(it),it) = SpolysIndep(1:Nindeppolys(it+1),it+1)
830	NpolysIndep(1:Nindeppolys(it),it) = NpolysIndep(1:Nindeppolys(it+1),it+1)
831	Npolystime(it) = Npolystime(it+1)
832
833	DO ip=1, Nindeppolys(it)
834	polysIndep(ip,1,it) = polysIndep(ip,1,it+1)
835	polysIndep(ip,2,it) = polysIndep(ip,2,it+1)
836	END DO
837	SpolysIndep(:,it+1) = 0
838	NpolysIndep(:,it+1) = 0
839	polysIndep(:,:,it+1) = 0
840
841	END DO
842	CLOSE(fprevunit)
843	IF (dbg) PRINT *," Succesful writting of ASCII chain of polygons 'polygons_overlap.dat' !!"
844	END IF
845	! ASCII file for all the polygons and their previous associated one
846	fprevunit = freeunit()
847	OPEN(fprevunit, file='polygons_overlap.dat', status='old', form='formatted', iostat=ios)
848	msg = "Problems opening file: 'polygons_overlap.dat'"
849	CALL ErrMsg(msg, fname, ios)
850
851	it = 1
852	IF (dbg) THEN
853	PRINT *, 'Coincidencies to connect _______'
854	DO iit=1, dt
855	! Reading from the ASCII file Independent polygons and their associated
856	CALL read_overlap_polys_ascii(fprevunit, iit, Nmaxpoly, Nindeppolys(it), SpolysIndep(:,it), &
857	NpolysIndep(:,it), polysIndep(:,:,it))
858	PRINT *,' it:', iit, ' Nindep:', Nindeppolys(it)
859	PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs'
860	DO ip=1, Nindeppolys(it)
861	PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':', &
862	polysIndep(ip,1:NpolysIndep(ip,it),it)
863	END DO
864	END DO
865	END IF
866
867	REWIND(fprevunit)
868
869	! Trajectories
870	! It should be done following the number of 'independent' polygons
871	! One would concatenate that independent polygons which share IDs from one step to another
872	IF (dotracks) THEN
873
874	! ASCII file for the trajectories
875	ftrackunit = freeunit()
876	OPEN(ftrackunit, file='trajectories_overlap.dat', status='new', form='formatted', iostat=ios)
877	msg = "Problems opening file: 'trajectories_overlap.dat'"
878	IF (ios == 17) PRINT *," Careful: 'trajectories_overlap.dat' already exists!!"
879	CALL ErrMsg(msg,fname,ios)
880
881	! First time-step. Take all polygons
882	itrack = 0
883	tracks = zeroRK
884	Ntracks = 0
885	it = 1
886	iit = 1
887	CALL read_overlap_polys_ascii(fprevunit, iit, Nmaxpoly, Nindeppolys(it), SpolysIndep(:,it), &
888	NpolysIndep(:,it), polysIndep(:,:,it))
889
890	DO ip=1, Nindeppolys(1)
891	itrack = itrack + 1
892	tracks(1,1,itrack,1) = itrack*1.
893	tracks(2,1,itrack,1) = SpolysIndep(ip,1)
894	tracks(3,1,itrack,1) = ctrpolys(1,ip,1)
895	tracks(4,1,itrack,1) = ctrpolys(2,ip,1)
896	tracks(5,1,itrack,1) = 1
897	Ntracks(1) = Ntracks(1) + 1
898	END DO
899
900	! Writting first time-step trajectories to the intermediate file
901	CALL write_overlap_tracks_ascii(ftrackunit,iit,Nmaxpoly, Ntracks(it), tracks(:,:,1:Ntracks(it),it))
902
903	! Looping allover already assigned tracks
904	it = 2
905	maxtrack = Ntracks(1)
906	timesteps: DO iit=2, dt
907	CALL read_overlap_polys_ascii(fprevunit, iit, Nmaxpoly, Nindeppolys(it), SpolysIndep(:,it), &
908	NpolysIndep(:,it), polysIndep(:,:,it))
909	IF (dbg) PRINT *,'track-timestep:', iit, 'N indep polys:', Nindeppolys(it)
910	! Indep polygons current time-step
911	current_poly: DO i=1, Nindeppolys(it)
912	IF (dbg) PRINT *,' curent poly:', i, 'Prev poly:', SpolysIndep(i,it), ' N ass. polygons:', &
913	NpolysIndep(i,it), 'ass. poly:', polysIndep(i,1:NpolysIndep(i,it),it)
914	Indeppolychained = .FALSE.
915
916	! Number of tracks previous time-step
917	! Looping overall
918	it1_tracks: DO itt=1, Ntracks(it-1)
919	itrack = tracks(1,1,itt,it-1)
920	! Number polygons ID assigned
921	Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0)
922	IF (dbg) PRINT *,itt,' track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1)
923
924	! Looking for coincidencies
925	DO iip=1, Ntprev
926	IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN
927	Indeppolychained = .TRUE.
928	IF (dbg) PRINT *,' coincidence found by polygon:', tracks(2,iip,itt,it-1)
929	EXIT
930	END IF
931	END DO
932	IF (Indeppolychained) THEN
933	Ntracks(it) = Ntracks(it) + 1
934	ictrack = Ntracks(it)
935	! Assigning all the IDs to the next step of the track
936	DO iip=1, NpolysIndep(i,it)
937	iiprev = polysIndep(i,iip,it)
938	tracks(1,iip,ictrack,it) = itrack
939	tracks(2,iip,ictrack,it) = iiprev
940	tracks(3,iip,ictrack,it) = ctrpolys(1,iiprev,iit)
941	tracks(4,iip,ictrack,it) = ctrpolys(2,iiprev,iit)
942	tracks(5,iip,ictrack,it) = iit
943	END DO
944	EXIT
945	END IF
946	IF (Indeppolychained) EXIT
947	END DO it1_tracks
948
949	! Creation of a new track
950	IF (.NOT.Indeppolychained) THEN
951	Ntracks(it) = Ntracks(it) + 1
952	ictrack = Ntracks(it)
953	! ID of new track
954	maxtrack = maxtrack + 1
955	IF (dbg) PRINT *,' New track!', maxtrack
956
957	! Assigning all the IDs to the next step of the track
958	DO j=1, NpolysIndep(i,it)
959	iiprev = polysIndep(i,j,it)
960	tracks(1,j,ictrack,it) = maxtrack
961	tracks(2,j,ictrack,it) = iiprev
962	tracks(3,j,ictrack,it) = ctrpolys(1,iiprev,iit)
963	tracks(4,j,ictrack,it) = ctrpolys(2,iiprev,iit)
964	tracks(5,j,ictrack,it) = iit
965	END DO
966	END IF
967
968	END DO current_poly
969
970	IF (dbg) THEN
971	PRINT *,' At this time-step:', iit, ' N trajectories:', Ntracks(it)
972	DO i=1, Ntracks(it)
973	Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0)
974	PRINT *,i ,'ID tracks:', tracks(1,1,i,it), 'ID polygon:', tracks(2,1:Nprev,i,it)
975	END DO
976	END IF
977
978	CALL write_overlap_tracks_ascii(ftrackunit,iit,Nmaxpoly,Ntracks(it),tracks(:,:,1:Ntracks(it),it))
979	! Re-initializing for the next time-step
980	tracks(:,:,:,it-1) = zeroRK
981	Ntracks(it-1) = Ntracks(it)
982	tracks(:,:,1:Ntracks(it-1),it-1) = tracks(:,:,1:Ntracks(it),it)
983	Ntracks(it) = 0
984	tracks(:,:,:,it) = zeroRK
985
986	END DO timesteps
987	CLOSE(ftrackunit)
988	IF (dbg) PRINT *," Succesful writting of ASCII chain of polygons 'trajectories_overlap.dat' !!"
989	CLOSE(fprevunit)
990	END IF
991
992	! Summarizing trajectories
993	! When multiple polygons are available, the mean of their central positions determines the position
994
995	IF (doftracks) THEN
996	! ASCII file for the trajectories
997	ftrackunit = freeunit()
998	OPEN(ftrackunit, file='trajectories_overlap.dat', status='old', form='formatted', iostat=ios)
999	msg = "Problems opening file: 'trajectories_overlap.dat'"
1000	CALL ErrMsg(msg,fname,ios)
1001
1002	! ASCII file for the final trajectories
1003	ftrunit = freeunit()
1004	OPEN(ftrunit, file='trajectories.dat', status='new', form='formatted', iostat=ios)
1005	msg = "Problems opening file: 'trajectories.dat'"
1006	IF (ios == 17) PRINT *," Careful: 'trajectories.dat' already exists!!"
1007	CALL ErrMsg(msg,fname,ios)
1008
1009	finaltracks = zeroRK
1010
1011	DO itt=1, Nmaxtracks
1012	CALL read_overlap_single_track_ascii(ftrackunit, dt, Nmaxpoly, Nmaxtracks, itt, singletrack)
1013
1014	! It might reach the las trajectory
1015	IF (ALL(singletrack == zeroRK)) EXIT
1016
1017	itrack = INT(MAXVAL(singletrack(1,1,:)))
1018	IF (dbg) THEN
1019	PRINT *,' Trajectory:', itt, '_______', itrack
1020	DO it=1, dt
1021	IF (singletrack(2,1,it) /= zeroRK) THEN
1022	j = COUNT(singletrack(2,:,it) /= zeroRK)
1023	PRINT *,it,':',(singletrack(3,i,it),',',singletrack(4,i,it),' ; ',i=1,j)
1024	END IF
1025	END DO
1026	END IF
1027
1028	finaltracks = zeroRK
1029	finaltracks(1,:) = itrack*oneRK
1030	DO it =1, dt
1031	Nprev = COUNT(INT(singletrack(2,:,it)) /= zeroRK)
1032	IF (Nprev /= 0) THEN
1033	finaltracks(5,it) = it*oneRK
1034	IF (TRIM(methodmulti) == 'largest') THEN
1035	maxarea = -10.*oneRK
1036	DO ip=1, Nprev
1037	IF (areapolys(singletrack(2,ip,it),it) > maxarea) THEN
1038	maxarea = areapolys(singletrack(2,ip,it),it)
1039	i = ip
1040	END IF
1041	END DO
1042	IF (dbg) THEN
1043	PRINT *,' Determine the trajectory coordinates to the largest polygon:', i, &
1044	' area:', maxarea
1045	END IF
1046	finaltracks(2,it) = singletrack(2,i,it)*oneRK
1047	finaltracks(3,it) = singletrack(3,i,it)
1048	finaltracks(4,it) = singletrack(4,i,it)
1049	ELSE IF (TRIM(methodmulti) == 'closest') THEN
1050	IF (it > 1) THEN
1051	mindist = 10000000.*oneRK
1052	DO ip=1, Nprev
1053	dist = SQRT((singletrack(3,ip,it)-finaltracks(3,it-1))**2 + &
1054	(singletrack(4,ip,it)-finaltracks(4,it-1))**2 )
1055	IF (dist < mindist) THEN
1056	mindist = dist
1057	i = ip
1058	END IF
1059	END DO
1060	finaltracks(2,it) = singletrack(3,i,it)*oneRK
1061	finaltracks(3,it) = singletrack(3,i,it)
1062	finaltracks(4,it) = singletrack(4,i,it)
1063	IF (dbg) THEN
1064	PRINT *,' Determine the trajectory coordinates to the closest previous polygon:',i,&
1065	' distance:', mindist
1066	END IF
1067	ELSE
1068	maxarea = -10.*oneRK
1069	DO ip=1, Nprev
1070	IF (areapolys(singletrack(2,ip,it),it) > maxarea) THEN
1071	maxarea = areapolys(singletrack(2,ip,it),it)
1072	i = ip
1073	END IF
1074	END DO
1075	IF (dbg) THEN
1076	PRINT *, ' Determine the trajectory coordinates to the largest polygon:', i, &
1077	' area:', maxarea, ' at the first time-step then to the closest'
1078	END IF
1079	finaltracks(2,it) = i*oneRK
1080	finaltracks(3,it) = singletrack(3,i,it)
1081	finaltracks(4,it) = singletrack(4,i,it)
1082	END IF
1083	ELSE
1084	totArea = zeroRK
1085	finaltracks(2,it) = -oneRK
1086	finaltracks(3,it) = zeroRK
1087	finaltracks(4,it) = zeroRK
1088	DO ip=1, Nprev
1089	areai = areapolys(singletrack(2,ip,it),it)
1090	totArea = totArea + areai
1091	finaltracks(3,it) = finaltracks(3,it) + singletrack(3,ip,it)*areai
1092	finaltracks(4,it) = finaltracks(4,it) + singletrack(4,ip,it)*areai
1093	END DO
1094	finaltracks(3,it) = finaltracks(3,it)/totArea
1095	finaltracks(4,it) = finaltracks(4,it)/totArea
1096	IF (dbg) THEN
1097	PRINT *,' Determine the trajectory coordinates to the area-averaged polygon ' // &
1098	' total area:', totArea
1099	END IF
1100
1101	END IF
1102
1103	END IF
1104	END DO
1105	! Writting the final track into the ASCII file
1106	CALL write_finaltrack_ascii(ftrunit, dt, finaltracks)
1107
1108	END DO
1109	CLOSE(ftrackunit)
1110	IF (dbg) PRINT *," Succesful writting of ASCII trajectories 'trajectories.dat' !!"
1111	CLOSE(ftrunit)
1112	END IF
1113
1114	IF (ALLOCATED(coins)) DEALLOCATE(coins)
1115	IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts)
1116
1117	RETURN
1118
1119	END SUBROUTINE poly_overlap_tracks_area_ascii
1120
1121	SUBROUTINE poly_overlap_tracks_area(dbg, dx, dy, dt, minarea, inNallpolys, allpolys, ctrpolys, &
1122	areapolys, Nmaxpoly, Nmaxtracks, tracks, finaltracks)
1123	! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum
1124	! overlaping/coincidence filtrered by a minimal area
1125
1126	IMPLICIT NONE
1127
1128	LOGICAL, INTENT(in) :: dbg
1129	INTEGER, INTENT(in) :: dx, dy, dt, Nmaxpoly, Nmaxtracks
1130	INTEGER, DIMENSION(dt), INTENT(in) :: inNallpolys
1131	INTEGER, DIMENSION(dx,dy,dt), INTENT(in) :: allpolys
1132	REAL(r_k), INTENT(in) :: minarea
1133	REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in) :: ctrpolys
1134	REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in) :: areapolys
1135	REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,dt), &
1136	INTENT(out) :: tracks
1137	REAL(r_k), DIMENSION(4,Nmaxtracks,dt), INTENT(out) :: finaltracks
1138
1139	! Local
1140	INTEGER :: i, j, ip, it, iip, itt
1141	INTEGER :: ierr
1142	REAL(r_k), DIMENSION(Nmaxpoly) :: Aprevpolys, Acurrpolys
1143	REAL(r_k), DIMENSION(2,Nmaxpoly) :: Cprevpolys, Ccurrpolys
1144	INTEGER, DIMENSION(dt) :: Nallpolys
1145	INTEGER, DIMENSION(dx,dy) :: meetpolys, searchpolys
1146	INTEGER, DIMENSION(Nmaxpoly) :: coincidencies
1147	INTEGER, DIMENSION(Nmaxpoly) :: prevID, currID
1148	INTEGER, DIMENSION(:), ALLOCATABLE :: coins
1149	INTEGER, DIMENSION(:,:), ALLOCATABLE :: coinsNpts
1150	INTEGER :: Nmeet, Nsearch, Nindep
1151	INTEGER, DIMENSION(dt) :: Nindeppolys
1152	CHARACTER(len=5) :: NcoinS
1153	INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt) :: polysIndep
1154	INTEGER, DIMENSION(Nmaxpoly,dt) :: NpolysIndep
1155	INTEGER, DIMENSION(Nmaxpoly,dt) :: SpolysIndep
1156	INTEGER :: iindep, iiprev
1157	INTEGER :: Nprev, NNprev, Ntprev
1158	LOGICAL :: Indeppolychained
1159	INTEGER :: itrack, ictrack
1160	REAL(r_k) :: ixp, iyp
1161	INTEGER :: ttrack
1162	INTEGER, DIMENSION(dt) :: Ntracks
1163	INTEGER :: idtrack, maxtrack
1164
1165	!!!!!!! Variables
1166	! dx,dy,dt: space/time dimensions
1167	! Nallpolys: Vector with the number of polygons at each time-step
1168	! allpolys: Series of 2D field with the polygons
1169	! minarea: minimal area (in same units as areapolys) to perform the tracking
1170	! ctrpolys: center of the polygons
1171	! areapolys: area of the polygons
1172	! Nmaxpoly: Maximum possible number of polygons
1173	! Nmaxtracks: maximum number of tracks
1174	! tracks: series of consecutive polygons
1175	! trackperiod: period of the track in time-steps
1176
1177	fname = 'poly_overlap_tracks_area'
1178
1179	IF (dbg) PRINT *,TRIM(fname)
1180
1181	! Number of independent polygons by time step
1182	Nindeppolys = 0
1183	! Number of polygons attached to each independent polygons by time step
1184	NpolysIndep = 0
1185	! ID of searching polygon attached to each independent polygons by time step
1186	SpolysIndep = 0
1187	! ID of polygons attached to each independent polygons by time step
1188	polysIndep = 0
1189	! ID of polygons from previous time-step
1190	prevID = 0
1191	! ID of polygons from current time-step
1192	currID = 0
1193
1194	! First time-step all are independent polygons
1195	it = 1
1196	Nmeet = inNallpolys(it)
1197	Nindeppolys(it) = Nmeet
1198	ip = 0
1199	meetpolys = allpolys(:,:,it)
1200	DO i=1, Nmeet
1201	IF (areapolys(i,it) >= minarea) THEN
1202	ip = ip + 1
1203	SpolysIndep(ip,it) = i
1204	currID(ip) = i
1205	Acurrpolys(ip) = areapolys(i,it)
1206	Ccurrpolys(1,ip) = ctrpolys(1,i,it)
1207	Ccurrpolys(2,ip) = ctrpolys(2,i,it)
1208	NpolysIndep(ip,it) = 1
1209	polysIndep(ip,1,it) = i
1210	ELSE
1211	WHERE (meetpolys == i)
1212	meetpolys = 0
1213	END WHERE
1214	END IF
1215	END DO
1216	Nallpolys(1) = ip
1217	Nindeppolys(1) = ip
1218
1219	! Starting step
1220	it = 0
1221	IF (dbg) THEN
1222	PRINT *,' time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',0
1223	PRINT *,' number of independent polygons:', Nindeppolys(it+1)
1224	PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______'
1225	DO i=1, Nindeppolys(it+1)
1226	PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', &
1227	polysIndep(i,1:NpolysIndep(i,it+1),it+1)
1228	END DO
1229	END IF
1230
1231	! Looking for the coincidencies at each time step
1232	DO it=1, dt-1
1233	! Number of times that a polygon has a coincidence
1234	coincidencies = 0
1235
1236	Nmeet = inNallpolys(it+1)
1237	searchpolys = meetpolys
1238	meetpolys = allpolys(:,:,it+1)
1239	prevID = 0
1240	prevID = currID
1241	Aprevpolys = Acurrpolys
1242	Cprevpolys = Ccurrpolys
1243	ip = 0
1244
1245	DO i=1, Nmeet
1246	IF (areapolys(i,it+1) >= minarea) THEN
1247	ip = ip + 1
1248	currID(ip) = i
1249	Acurrpolys(ip) = areapolys(i,it+1)
1250	Acurrpolys(ip) = areapolys(i,it+1)
1251	Ccurrpolys(1,ip) = ctrpolys(1,i,it+1)
1252	Ccurrpolys(2,ip) = ctrpolys(2,i,it+1)
1253	ELSE
1254	WHERE (meetpolys == i)
1255	meetpolys = 0
1256	END WHERE
1257	END IF
1258	END DO
1259	Nallpolys(it+1) = ip
1260	Nindeppolys(it+1) = ip
1261
1262	Nmeet = Nallpolys(it+1)
1263	! Looking throughout the independent polygons
1264	Nsearch = Nindeppolys(it)
1265
1266	IF (ALLOCATED(coins)) DEALLOCATE(coins)
1267	ALLOCATE(coins(Nmeet), STAT=ierr)
1268	msg="Problems allocating 'coins'"
1269	CALL ErrMsg(msg,fname,ierr)
1270
1271	IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts)
1272	ALLOCATE(coinsNpts(Nmeet, Nsearch), STAT=ierr)
1273	msg="Problems allocating 'coinsNpts'"
1274	CALL ErrMsg(msg,fname,ierr)
1275
1276	CALL coincidence_all_polys_area(dbg, dx,dy, Nmeet, currID, meetpolys, Acurrpolys(1:Nmeet), &
1277	Nsearch, prevID, searchpolys, Cprevpolys(:,1:Nsearch), Aprevpolys(1:Nsearch), coins, &
1278	coinsNpts)
1279
1280	! Counting the number of times that a polygon has a coincidency
1281	IF (dbg) THEN
1282	PRINT *,' Coincidencies for the given time-step:', it+1,' _______'
1283	DO i=1, Nmeet
1284	PRINT *,currID(i), coins(i),' N search pts:', coinsNpts(i,:)
1285	END DO
1286	END IF
1287
1288	! Looking for the same equivalencies
1289	Nindep = 0
1290	DO i=1, Nmeet
1291	IF (coins(i) == -1) THEN
1292	Nindep = Nindep + 1
1293	SpolysIndep(Nindep,it+1) = -1
1294	NpolysIndep(Nindep,it+1) = NpolysIndep(Nindep,it+1) + 1
1295	polysIndep(Nindep,NpolysIndep(Nindep,it+1),it+1) = currID(i)
1296	ELSE IF (coins(i) == -9) THEN
1297	WRITE(NcoinS,'(I5)')coins(i)
1298	msg="coins= "//TRIM(NcoinS)//" This is an error. One should have always only one " // &
1299	"coincidence of polygon"
1300	CALL ErrMsg(msg, fname, -1)
1301	ELSE
1302	! Looking for coincidences with previous independent polygons
1303	DO ip=1, Nsearch
1304	! Looking into the polygons associated
1305	NNprev = NpolysIndep(ip,it)
1306	DO j=1, NNprev
1307	IF (coins(i) == polysIndep(ip,j,it)) THEN
1308	! Which index corresponds to this coincidence?
1309	iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
1310	IF (iindep == -1) THEN
1311	Nindep = Nindep + 1
1312	SpolysIndep(Nindep,it+1) = coins(i)
1313	END IF
1314	iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
1315	IF (iindep < 0) THEN
1316	PRINT *,' Looking for:', coins(i)
1317	PRINT *,' Within:', SpolysIndep(1:Nindep,it+1)
1318	PRINT *,' Might content:', polysIndep(ip,1:NNprev,it)
1319	PRINT *,' From an initial list:', coins(1:Nmeet)
1320	msg = 'Wrong index! There must be an index here'
1321	CALL ErrMsg(msg,fname,iindep)
1322	END IF
1323	coincidencies(ip) = coincidencies(ip) + 1
1324	NpolysIndep(iindep,it+1) = NpolysIndep(iindep,it+1) + 1
1325	polysIndep(iindep,NpolysIndep(iindep,it+1),it+1) = currID(i)
1326	EXIT
1327	END IF
1328	END DO
1329	END DO
1330	END IF
1331	END DO
1332	Nindeppolys(it+1) = Nindep
1333
1334	IF (dbg) THEN
1335	PRINT *,' time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',Nsearch
1336	PRINT *,' number of independent polygons:', Nindeppolys(it+1)
1337	PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______'
1338	DO i=1, Nindeppolys(it+1)
1339	PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', &
1340	polysIndep(i,1:NpolysIndep(i,it+1),it+1)
1341	END DO
1342	END IF
1343	END DO
1344
1345	IF (dbg) THEN
1346	PRINT *, 'Coincidencies to connect _______'
1347	DO it=1, dt
1348	PRINT *,' it:', it, ' Nindep:', Nindeppolys(it)
1349	PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs'
1350	DO ip=1, Nindeppolys(it)
1351	PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':', &
1352	polysIndep(ip,1:NpolysIndep(ip,it),it)
1353	END DO
1354	END DO
1355
1356	END IF
1357
1358	! Trajectories
1359	! It should be done following the number of 'independent' polygons
1360	! One would concatenate that independent polygons which share IDs from one step to another
1361
1362	! First time-step. Take all polygons
1363	itrack = 0
1364	tracks = 0.
1365	Ntracks = 0
1366	DO ip=1, Nindeppolys(1)
1367	itrack = itrack + 1
1368	tracks(1,1,itrack,1) = itrack*1.
1369	tracks(2,1,itrack,1) = SpolysIndep(ip,1)
1370	tracks(3,1,itrack,1) = ctrpolys(1,ip,1)
1371	tracks(4,1,itrack,1) = ctrpolys(2,ip,1)
1372	tracks(5,1,itrack,1) = 1
1373	Ntracks(1) = Ntracks(1) + 1
1374	END DO
1375
1376	! Looping allover already assigned tracks
1377	timesteps: DO it=2, dt
1378	IF (dbg) PRINT *,'track-timestep:', it, 'N indep polys:', Nindeppolys(it)
1379	! Indep polygons current time-step
1380	current_poly: DO i=1, Nindeppolys(it)
1381	IF (dbg) PRINT *,' curent poly:', i, 'Prev poly:', SpolysIndep(i,it), ' N ass. polygons:', &
1382	NpolysIndep(i,it), 'ass. poly:', polysIndep(i,1:NpolysIndep(i,it),it)
1383	Indeppolychained = .FALSE.
1384
1385	! Number of tracks previous time-step
1386	! Looping overall
1387	it1_tracks: DO itt=1, Ntracks(it-1)
1388	itrack = tracks(1,1,itt,it-1)
1389	! Number polygons ID assigned
1390	Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0)
1391	IF (dbg) PRINT *,itt,' track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1)
1392
1393	! Looking for coincidencies
1394	DO iip=1, Ntprev
1395	IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN
1396	Indeppolychained = .TRUE.
1397	IF (dbg) PRINT *,' coincidence found by polygon:', tracks(2,iip,itt,it-1)
1398	EXIT
1399	END IF
1400	END DO
1401	IF (Indeppolychained) THEN
1402	Ntracks(it) = Ntracks(it) + 1
1403	ictrack = Ntracks(it)
1404	! Assigning all the IDs to the next step of the track
1405	DO iip=1, NpolysIndep(i,it)
1406	iiprev = polysIndep(i,iip,it)
1407	tracks(1,iip,ictrack,it) = itrack
1408	tracks(2,iip,ictrack,it) = iiprev
1409	ixp = ctrpolys(1,iiprev,it)
1410	iyp = ctrpolys(2,iiprev,it)
1411	tracks(3,iip,ictrack,it) = ixp
1412	tracks(4,iip,ictrack,it) = iyp
1413	tracks(5,iip,ictrack,it) = it
1414	END DO
1415	EXIT
1416	END IF
1417	IF (Indeppolychained) EXIT
1418	END DO it1_tracks
1419
1420	! Creation of a new track
1421	IF (.NOT.Indeppolychained) THEN
1422	Ntracks(it) = Ntracks(it) + 1
1423	ictrack = Ntracks(it)
1424	! ID of new track
1425	maxtrack = INT(MAXVAL(tracks(1,:,:,:)*1.))
1426	IF (dbg) PRINT *,' New track!', maxtrack+1
1427
1428	! Assigning all the IDs to the next step of the track
1429	DO j=1, NpolysIndep(i,it)
1430	iiprev = polysIndep(i,j,it)
1431	tracks(1,j,ictrack,it) = maxtrack+1
1432	tracks(2,j,ictrack,it) = iiprev
1433	ixp = ctrpolys(1,iiprev,it)
1434	iyp = ctrpolys(2,iiprev,it)
1435	tracks(3,j,ictrack,it) = ixp
1436	tracks(4,j,ictrack,it) = iyp
1437	tracks(5,j,ictrack,it) = it
1438	END DO
1439	END IF
1440
1441	END DO current_poly
1442
1443	IF (dbg) THEN
1444	PRINT *,' At this time-step:', it, ' N trajectories:', Ntracks(it)
1445	DO i=1, Ntracks(it)
1446	Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0)
1447	PRINT *,i ,'ID tracks:', tracks(1,1,i,it), 'ID polygon:', tracks(2,1:Nprev,i,it)
1448	END DO
1449	END IF
1450
1451	END DO timesteps
1452
1453	! Summarizing trajectories
1454	! When multiple polygons are available, the mean of their central positions determines the position
1455
1456	finaltracks = 0.
1457	maxtrack = MAXVAL(tracks(1,:,:,:))
1458
1459	DO it=1, dt
1460	DO itt=1, Ntracks(it)
1461	itrack = INT(tracks(1,1,itt,it))
1462	Nprev = COUNT(INT(tracks(2,:,itt,it)) /= 0)
1463	finaltracks(1,itrack,it) = itrack*1.
1464	finaltracks(2,itrack,it) = SUM(tracks(3,:,itt,it))/Nprev*1.
1465	finaltracks(3,itrack,it) = SUM(tracks(4,:,itt,it))/Nprev*1.
1466	finaltracks(4,itrack,it) = it*1.
1467	END DO
1468	END DO
1469
1470	DEALLOCATE(coins)
1471	DEALLOCATE(coinsNpts)
1472
1473	RETURN
1474
1475	END SUBROUTINE poly_overlap_tracks_area
1476
1477	SUBROUTINE coincidence_all_polys_area(dbg, dx, dy, Nallpoly, IDallpoly, allpoly, icpolys, Npoly, &
1478	IDpolys, polys, cpolys, apolys, polycoins, coinNptss)
1479	! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
1480	! In case of multiple coincidencies, the closest and then the largest is taken filtrered by a minimal area
1481	! Here the difference is that the index does not coincide with its ID
1482
1483	IMPLICIT NONE
1484
1485	LOGICAL, INTENT(in) :: dbg
1486	INTEGER, INTENT(in) :: dx, dy, Nallpoly, Npoly
1487	INTEGER, DIMENSION(dx,dy), INTENT(in) :: allpoly, polys
1488	INTEGER, DIMENSION(Nallpoly), INTENT(in) :: IDallpoly
1489	INTEGER, DIMENSION(Npoly), INTENT(in) :: IDpolys
1490	REAL(r_k), DIMENSION(2,Nallpoly), INTENT(in) :: icpolys
1491	REAL(r_k), DIMENSION(2,Npoly), INTENT(in) :: cpolys
1492	REAL(r_k), DIMENSION(Npoly), INTENT(in) :: apolys
1493	INTEGER, DIMENSION(Nallpoly), INTENT(out) :: polycoins
1494	INTEGER, DIMENSION(Nallpoly,Npoly), INTENT(out) :: coinNptss
1495
1496	! Local
1497	INTEGER :: i, j, ip
1498	INTEGER :: maxcorr
1499	INTEGER :: Nmaxcorr
1500	LOGICAL, DIMENSION(dx,dy) :: boolpoly
1501	INTEGER :: maxcoin
1502	REAL :: dist, maxcoindist, maxcoinarea
1503	INTEGER, DIMENSION(Npoly) :: IDcoins
1504
1505	!!!!!!! Variables
1506	! dx,dy: dimension of the space
1507	! Nallpoly: Number of polygons to find coincidence
1508	! allpoly: space with the polygons to meet
1509	! IDallpoly: ID of the polygon to find coincidence
1510	! icpolys: center of the polygons to look for the coincidence
1511	! Npoly: number of polygons on the 2D space
1512	! polys: 2D field of polygons identified by their integer number (0 for no polygon)
1513	! IDpolys: ID of the polygon to search for coincidences
1514	! cpolys: center of the polygons
1515	! apolys: area of the polygons
1516	! polycoins: coincident polyogn
1517	! -1: no-coincidence
1518	! 1 < Npoly: single coincidence with a given polygon
1519	! -9: coincidence with more than one polygon
1520	! coinNptss: number of points coincident with each polygon
1521
1522	fname = 'coincidence_all_polys_area'
1523	IF (dbg) PRINT *,TRIM(fname)
1524
1525	DO ip=1, Nallpoly
1526	boolpoly = allpoly == IDallpoly(ip)
1527	CALL coincidence_poly_area(dbg, dx, dy, boolpoly, Npoly, polys, polycoins(ip), IDcoins, &
1528	coinNptss(ip,:))
1529	IF (dbg) PRINT *,' polygon', IDallpoly(ip), ' coincidence with:', polycoins(ip), 'IDpolys:', IDpolys(1:Npoly)
1530
1531	! Coincidence with more than one polygon
1532	IF (polycoins(ip) == -9) THEN
1533	maxcoindist = -10.
1534	maxcoinarea = -10.
1535	maxcoin = MAXVAL(coinNptss(ip,:))
1536	DO j=1, Npoly
1537	IF (coinNptss(ip,j) == maxcoin) THEN
1538	dist = SQRT( (icpolys(1,ip)1.-cpolys(1,j)1.)*2 + (icpolys(2,ip)1.-cpolys(2,j)1.)*2 )
1539	IF ( dist > maxcoindist) THEN
1540	maxcoindist = dist
1541	maxcoinarea = apolys(j)
1542	polycoins(ip) = IDcoins(j)
1543	ELSE IF ( dist == maxcoindist) THEN
1544	IF (apolys(j) > maxcoinarea) THEN
1545	polycoins(ip) = IDcoins(j)
1546	maxcoinarea = apolys(j)
1547	END IF
1548	END IF
1549	END IF
1550	END DO
1551	END IF
1552	END DO
1553
1554	RETURN
1555
1556	END SUBROUTINE coincidence_all_polys_area
1557
1558	SUBROUTINE coincidence_poly_area(dbg, dx, dy, poly, Npoly, polys, polycoin, IDpoly, coinNpts)
1559	! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
1560	! Here the difference is that the index does not coincide with its ID
1561
1562	IMPLICIT NONE
1563
1564	LOGICAL, INTENT(in) :: dbg
1565	INTEGER, INTENT(in) :: dx, dy, Npoly
1566	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: poly
1567	INTEGER, DIMENSION(dx,dy), INTENT(in) :: polys
1568	INTEGER, INTENT(out) :: polycoin
1569	INTEGER, DIMENSION(Npoly), INTENT(out) :: IDpoly, coinNpts
1570
1571	! Local
1572	INTEGER :: i, j, ip
1573	INTEGER :: maxcorr
1574	INTEGER :: Nmaxcorr
1575	! Lluis
1576	INTEGER :: Ndiffvals
1577	INTEGER, DIMENSION(:), ALLOCATABLE :: diffvals
1578
1579	!!!!!!! Variables
1580	! dx,dy: dimension of the space
1581	! poly: bolean polygon to meet
1582	! Npoly: number of polygons on the 2D space
1583	! polys: 2D field of polygons identified by their integer number (0 for no polygon)
1584	! polycoin: coincident polyogn
1585	! -1: no-coincidence
1586	! 1 < Npoly: single coincidence with a given polygon
1587	! -9: coincidence with more than one polygon
1588	! IDpoly: ID of the found polygon
1589	! coinNpts: number of points coincident with each polygon
1590
1591	fname = 'coincidence_poly_area'
1592	IF (dbg) PRINT *,TRIM(fname)
1593
1594	IF (dbg) THEN
1595	PRINT *,' Boolean polygon to search coincidences ...'
1596	DO i=1,dx
1597	PRINT *,poly(i,:)
1598	END DO
1599
1600	PRINT *,' 2D polygons space ...'
1601	DO i=1,dx
1602	PRINT '(1000(I7,1x))',polys(i,:)
1603	END DO
1604	END IF
1605
1606	IF (ALLOCATED(diffvals)) DEALLOCATE(diffvals)
1607	ALLOCATE(diffvals(dx*dy))
1608
1609	! Checking for consistency on number of polygons and real content (except 0 value)
1610	CALL Nvalues_2DArrayI(dx, dy, dx*dy, polys, Ndiffvals, diffvals)
1611	IF (Ndiffvals -1 /= Npoly) THEN
1612	PRINT *,TRIM(emsg)
1613	PRINT *,' number of different values:', Ndiffvals-1, ' theoretical Npoly:', Npoly
1614	PRINT *,' Different values:', diffvals(1:Ndiffvals)
1615	msg = 'Number of different values and Npoly must coincide'
1616	CALL ErrMsg(msg, fname, -1)
1617	END IF
1618
1619	! Looking for coincient points for the polygon
1620	coinNpts = 0
1621	IDpoly = 0
1622	ip = 0
1623	DO i=1,dx
1624	DO j=1,dy
1625	IF (poly(i,j) .AND. polys(i,j) .NE. 0) THEN
1626	IF (.NOT.ANY(IDpoly == polys(i,j))) THEN
1627	ip = ip + 1
1628	IDpoly(ip) = polys(i,j)
1629	ELSE
1630	ip = Index1DarrayI(IDpoly, Npoly, polys(i,j))
1631	END IF
1632	coinNpts(ip) = coinNpts(ip) + 1
1633	END IF
1634	END DO
1635	END DO
1636
1637	maxcorr = 0
1638	maxcorr = INT(MAXVAL(coinNpts*1.))
1639
1640	IF (dbg) PRINT *,' Maximum coincidence:', maxcorr
1641	IF (maxcorr == 0) THEN
1642	polycoin = -1
1643	ELSE
1644	Nmaxcorr = 0
1645	DO ip=1, Npoly
1646	IF (coinNpts(ip) == maxcorr) THEN
1647	Nmaxcorr = Nmaxcorr+1
1648	polycoin = IDpoly(ip)
1649	END IF
1650	END DO
1651	IF (Nmaxcorr > 1) polycoin = -9
1652	END IF
1653
1654	IF (dbg) THEN
1655	PRINT *,' Coincidences for each polygon _______', Npoly
1656	DO ip=1, Npoly
1657	PRINT *,' ',ip, ' ID:', IDpoly(ip),': ', coinNpts(ip)
1658	END DO
1659	END IF
1660
1661	RETURN
1662
1663	END SUBROUTINE coincidence_poly_area
1664
1665	SUBROUTINE poly_overlap_tracks(dbg, dx, dy, dt, minarea, Nallpolys, allpolys, ctrpolys, &
1666	areapolys, Nmaxpoly, Nmaxtracks, tracks, finaltracks)
1667	! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum overlaping/coincidence
1668
1669	IMPLICIT NONE
1670
1671	LOGICAL, INTENT(in) :: dbg
1672	INTEGER, INTENT(in) :: dx, dy, dt, Nmaxpoly, Nmaxtracks
1673	INTEGER, DIMENSION(dt), INTENT(in) :: Nallpolys
1674	INTEGER, DIMENSION(dx,dy,dt), INTENT(in) :: allpolys
1675	REAL(r_k), INTENT(in) :: minarea
1676	REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in) :: ctrpolys
1677	REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in) :: areapolys
1678	REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,dt), &
1679	INTENT(out) :: tracks
1680	REAL(r_k), DIMENSION(4,Nmaxtracks,dt), INTENT(out) :: finaltracks
1681
1682	! Local
1683	INTEGER :: i, j, ip, it, iip, itt
1684	INTEGER :: ierr
1685	INTEGER, DIMENSION(Nmaxpoly,dt) :: coincidencies, NOcoincidencies
1686	INTEGER, DIMENSION(:), ALLOCATABLE :: coins
1687	INTEGER, DIMENSION(:,:), ALLOCATABLE :: coinsNpts
1688	INTEGER, DIMENSION(Nmaxpoly,dt) :: polycoincidencies
1689	INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt) :: coincidenciesNpts
1690	INTEGER :: Nmeet, Nsearch, Nindep
1691	INTEGER, DIMENSION(dt) :: Nindeppolys
1692	CHARACTER(len=5) :: NcoinS
1693	INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt) :: polysIndep
1694	INTEGER, DIMENSION(Nmaxpoly,dt) :: NpolysIndep
1695	INTEGER, DIMENSION(Nmaxpoly,dt) :: SpolysIndep
1696	INTEGER :: iindep, iiprev
1697	INTEGER :: Nprev, NNprev, Ntprev
1698	LOGICAL :: Indeppolychained
1699	INTEGER :: itrack, ictrack
1700	INTEGER :: ixp, iyp, ttrack
1701	INTEGER, DIMENSION(dt) :: Ntracks
1702	INTEGER :: idtrack, maxtrack
1703
1704	!!!!!!! Variables
1705	! dx,dy,dt: space/time dimensions
1706	! Nallpolys: Vector with the number of polygons at each time-step
1707	! allpolys: Series of 2D field with the polygons
1708	! minarea: minimal area (in same units as areapolys) to perform the tracking
1709	! ctrpolys: center of the polygons
1710	! areapolys: area of the polygons
1711	! Nmaxpoly: Maximum possible number of polygons
1712	! Nmaxtracks: maximum number of tracks
1713	! tracks: series of consecutive polygons
1714	! trackperiod: period of the track in time-steps
1715
1716	fname = 'poly_overlap_tracks'
1717
1718	IF (dbg) PRINT *,TRIM(fname)
1719
1720	polycoincidencies = fillvalI
1721	coincidenciesNpts = fillvalI
1722	! Number of times that a polygon has a coincidence
1723	coincidencies = 0
1724	! Polygons without a coincidence
1725	NOcoincidencies = 0
1726	! Number of independent polygons by time step
1727	Nindeppolys = 0
1728	! Number of polygons attached to each independent polygons by time step
1729	NpolysIndep = 0
1730	! ID of searching polygon attached to each independent polygons by time step
1731	SpolysIndep = 0
1732	! ID of polygons attached to each independent polygons by time step
1733	polysIndep = 0
1734
1735	! First time-step all are independent polygons
1736	it = 1
1737	Nmeet = Nallpolys(it)
1738	Nindeppolys(it) = Nmeet
1739	DO i=1, Nmeet
1740	SpolysIndep(i,it) = i
1741	NpolysIndep(1:Nmeet,it) = 1
1742	polysIndep(1,i,it) = i
1743	END DO
1744
1745	! Looking for the coincidencies at each time step
1746	DO it=1, dt-1
1747	Nmeet = Nallpolys(it+1)
1748	Nsearch = Nallpolys(it)
1749
1750	IF (ALLOCATED(coins)) DEALLOCATE(coins)
1751	ALLOCATE(coins(Nmeet), STAT=ierr)
1752	msg="Problems allocating 'coins'"
1753	CALL ErrMsg(msg,fname,ierr)
1754
1755	IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts)
1756	ALLOCATE(coinsNpts(Nmeet, Nsearch), STAT=ierr)
1757	msg="Problems allocating 'coinsNpts'"
1758	CALL ErrMsg(msg,fname,ierr)
1759
1760	CALL coincidence_all_polys(dbg, dx, dy, Nmeet, allpolys(:,:,it+1), ctrpolys(:,1:Nmeet,it+1), &
1761	Nsearch, allpolys(:,:,it), ctrpolys(:,1:Nsearch,it), areapolys(1:Nsearch,it), coins, coinsNpts)
1762
1763	polycoincidencies(1:Nmeet,it+1) = coins
1764	coincidenciesNpts(1:Nmeet,1:Nsearch,it+1) = coinsNpts
1765
1766	! Counting the number of times that a polygon has a coincidency
1767	IF (dbg) THEN
1768	PRINT *,' Coincidencies for the given time-step:', it+1,' _______'
1769	DO i=1, Nmeet
1770	PRINT *,coins(i),' N search pts:', coinsNpts(i,:)
1771	END DO
1772	END IF
1773
1774	Nindep = 0
1775	DO i=1, Nmeet
1776	IF (coins(i) == -1) THEN
1777	Nindep = Nindep + 1
1778	NOcoincidencies(i,it+1) = 1
1779	SpolysIndep(Nindep,it+1) = -1
1780	NpolysIndep(Nindep,it+1) = NpolysIndep(Nindep,it+1) + 1
1781	polysIndep(Nindep,NpolysIndep(Nindep,it+1),it+1) = i
1782	ELSE IF (coins(i) == -9) THEN
1783	WRITE(NcoinS,'(I5)')coins(i)
1784	msg="coins= "//TRIM(NcoinS)//" This is an error. One should have always only one " // &
1785	"coincidence of polygon"
1786	CALL ErrMsg(msg, fname, -1)
1787	ELSE
1788	DO ip=1, Nsearch
1789	IF (coins(i) == ip) THEN
1790	IF (coincidencies(ip,it+1) == 0) THEN
1791	Nindep = Nindep + 1
1792	SpolysIndep(Nindep,it+1) = ip
1793	END IF
1794	coincidencies(ip,it+1) = coincidencies(ip,it+1) + 1
1795	DO iindep=1, Nindep
1796	IF (SpolysIndep(iindep,it+1) == coins(i)) THEN
1797	NpolysIndep(iindep,it+1) = NpolysIndep(iindep,it+1) + 1
1798	polysIndep(iindep,NpolysIndep(iindep,it+1),it+1) = i
1799	END IF
1800	END DO
1801	END IF
1802	END DO
1803	END IF
1804	END DO
1805	Nindeppolys(it+1) = Nindep
1806
1807	IF (dbg) THEN
1808	PRINT *,' time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',Nsearch
1809	PRINT *,' number of independent polygons:', Nindeppolys(it+1)
1810	PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______'
1811	DO i=1, Nindeppolys(it+1)
1812	PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', &
1813	polysIndep(i,1:NpolysIndep(i,it+1),it+1)
1814	END DO
1815	END IF
1816	END DO
1817
1818	IF (dbg) THEN
1819	PRINT *, 'Coincidencies to connect _______'
1820	DO it=1, dt
1821	PRINT *,' it:', it, ' Nindep:', Nindeppolys(it)
1822	PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs'
1823	DO ip=1, Nindeppolys(it)
1824	PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':', &
1825	polysIndep(ip,1:NpolysIndep(ip,it),it)
1826	END DO
1827	END DO
1828
1829	END IF
1830
1831	! Trajectories
1832	! It should be done following the number of 'independent' polygons
1833	! One would concatenate that independent polygons which share IDs from one step to another
1834
1835	! First time-step. Take all polygons
1836	itrack = 0
1837	tracks = 0.
1838	Ntracks = 0
1839	DO ip=1, Nindeppolys(1)
1840	itrack = itrack + 1
1841	tracks(1,1,itrack,1) = itrack*1.
1842	tracks(2,1,itrack,1) = SpolysIndep(ip,1)
1843	tracks(3,1,itrack,1) = ctrpolys(1,ip,1)
1844	tracks(4,1,itrack,1) = ctrpolys(2,ip,1)
1845	tracks(5,1,itrack,1) = 1
1846	Ntracks(1) = Ntracks(1) + 1
1847	END DO
1848
1849	! Looping allover already assigned tracks
1850	timesteps: DO it=2, dt
1851	IF (dbg) PRINT *,'timestep:', it, 'N indep polys:', Nindeppolys(it)
1852	! Indep polygons current time-step
1853	current_poly: DO i=1, Nindeppolys(it)
1854	IF (dbg) PRINT *,' curent poly:', i, 'Prev poly:', SpolysIndep(i,it), ' N ass. polygons:', &
1855	NpolysIndep(i,it), 'ass. poly:', polysIndep(i,1:NpolysIndep(i,it),it)
1856	Indeppolychained = .FALSE.
1857
1858	! Number of tracks previous time-step
1859	! Looping overall
1860	it1_tracks: DO itt=1, Ntracks(it-1)
1861	itrack = tracks(1,1,itt,it-1)
1862	! Number polygons ID assigned
1863	Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0)
1864	IF (dbg) PRINT *,itt,' track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1)
1865
1866	! Looking for coincidencies
1867	DO iip=1, Ntprev
1868	IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN
1869	Indeppolychained = .TRUE.
1870	IF (dbg) PRINT *,' coincidence found by polygon:', tracks(2,iip,itt,it-1)
1871	EXIT
1872	END IF
1873	END DO
1874	IF (Indeppolychained) THEN
1875	Ntracks(it) = Ntracks(it) + 1
1876	ictrack = Ntracks(it)
1877	! Assigning all the IDs to the next step of the track
1878	DO iip=1, NpolysIndep(i,it)
1879	iiprev = polysIndep(i,iip,it)
1880	tracks(1,iip,ictrack,it) = itrack
1881	tracks(2,iip,ictrack,it) = iiprev
1882	ixp = ctrpolys(1,iiprev,it)
1883	iyp = ctrpolys(2,iiprev,it)
1884	tracks(3,iip,ictrack,it) = ixp
1885	tracks(4,iip,ictrack,it) = iyp
1886	tracks(5,iip,ictrack,it) = it
1887	END DO
1888	EXIT
1889	END IF
1890	END DO it1_tracks
1891
1892	! Creation of a new track
1893	IF (.NOT.Indeppolychained) THEN
1894	Ntracks(it) = Ntracks(it) + 1
1895	ictrack = Ntracks(it)
1896	! ID of new track
1897	maxtrack = INT(MAXVAL(tracks(1,:,:,:)*1.))
1898	IF (dbg) PRINT *,' New track!', maxtrack+1
1899
1900	! Assigning all the IDs to the next step of the track
1901	DO j=1, NpolysIndep(i,it)
1902	iiprev = polysIndep(i,j,it)
1903	tracks(1,j,ictrack,it) = maxtrack+1
1904	tracks(2,j,ictrack,it) = iiprev
1905	ixp = ctrpolys(1,iiprev,it)
1906	iyp = ctrpolys(2,iiprev,it)
1907	tracks(3,j,ictrack,it) = ixp
1908	tracks(4,j,ictrack,it) = iyp
1909	tracks(5,j,ictrack,it) = it
1910	END DO
1911	END IF
1912
1913	END DO current_poly
1914
1915	IF (dbg) THEN
1916	PRINT *,' At this time-step:', it, ' N trajectories:', Ntracks(it)
1917	DO i=1, Ntracks(it)
1918	Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0)
1919	PRINT *,i,'tracks:', tracks(2,1:Nprev,i,it)
1920	END DO
1921	END IF
1922
1923	END DO timesteps
1924
1925	! Summarizing trajectories
1926	! When multiple polygons are available, the mean of their central positions determines the position
1927
1928	finaltracks = 0.
1929	maxtrack = MAXVAL(tracks(1,:,:,:))
1930
1931	DO it=1, dt
1932	DO itt=1, Ntracks(it)
1933	itrack = INT(tracks(1,1,itt,it))
1934	Nprev = COUNT(INT(tracks(2,:,itt,it)) /= 0)
1935	PRINT *,'it:', it,'itrack:', itrack, 'Nprev:', Nprev
1936	finaltracks(1,itrack,it) = itrack*1.
1937	finaltracks(2,itrack,it) = SUM(tracks(3,:,itt,it))/Nprev
1938	finaltracks(3,itrack,it) = SUM(tracks(4,:,itt,it))/Nprev
1939	finaltracks(4,itrack,it) = it*1.
1940	PRINT *,' finaltrack:', finaltracks(:,itrack,it)
1941	END DO
1942	END DO
1943
1944	RETURN
1945
1946	END SUBROUTINE poly_overlap_tracks
1947
1948	SUBROUTINE coincidence_all_polys(dbg, dx, dy, Nallpoly, allpoly, icpolys, Npoly, polys, cpolys, &
1949	apolys, polycoins, coinNptss)
1950	! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
1951	! In case of multiple coincidencies, the closest and then the largest is taken
1952
1953	IMPLICIT NONE
1954
1955	LOGICAL, INTENT(in) :: dbg
1956	INTEGER, INTENT(in) :: dx, dy, Nallpoly, Npoly
1957	INTEGER, DIMENSION(dx,dy), INTENT(in) :: allpoly, polys
1958	REAL(r_k), DIMENSION(2,Nallpoly), INTENT(in) :: icpolys
1959	REAL(r_k), DIMENSION(2,Npoly), INTENT(in) :: cpolys
1960	REAL(r_k), DIMENSION(Npoly), INTENT(in) :: apolys
1961	INTEGER, DIMENSION(Nallpoly), INTENT(out) :: polycoins
1962	INTEGER, DIMENSION(Nallpoly,Npoly), INTENT(out) :: coinNptss
1963
1964	! Local
1965	INTEGER :: i, j, ip
1966	INTEGER :: maxcorr
1967	INTEGER :: Nmaxcorr
1968	LOGICAL, DIMENSION(dx,dy) :: boolpoly
1969	INTEGER :: maxcoin
1970	REAL :: dist, maxcoindist, maxcoinarea
1971
1972	!!!!!!! Variables
1973	! dx,dy: dimension of the space
1974	! Nallpoly: Number of polygons to find coincidence
1975	! allpoly: space with the polygons to meet
1976	! icpolys: center of the polygons to look for the coincidence
1977	! Npoly: number of polygons on the 2D space
1978	! polys: 2D field of polygons identified by their integer number (0 for no polygon)
1979	! cpolys: center of the polygons
1980	! apolys: area of the polygons
1981	! polycoins: coincident polyogn
1982	! -1: no-coincidence
1983	! 1 < Npoly: single coincidence with a given polygon
1984	! -9: coincidence with more than one polygon
1985	! coinNptss: number of points coincident with each polygon
1986
1987	fname = 'coincidence_all_polys'
1988	IF (dbg) PRINT *,TRIM(fname)
1989
1990	DO ip=1, Nallpoly
1991	boolpoly = allpoly == ip
1992	CALL coincidence_poly(dbg, dx, dy, boolpoly, Npoly, polys, polycoins(ip), coinNptss(ip,:))
1993	IF (dbg) PRINT *,' polygon', ip, ' coincidence with:', polycoins(ip)
1994
1995	! Coincidence with more than one polygon
1996	IF (polycoins(ip) == -9) THEN
1997	maxcoindist = -10.
1998	maxcoinarea = -10.
1999	maxcoin = MAXVAL(coinNptss(ip,:))
2000	DO j=1, Npoly
2001	IF (coinNptss(ip,j) == maxcoin) THEN
2002	dist = SQRT( (icpolys(1,ip)1.-cpolys(1,j)1.)*2 + (icpolys(2,ip)1.-cpolys(2,j)1.)*2 )
2003	IF ( dist > maxcoindist) THEN
2004	maxcoindist = dist
2005	maxcoinarea = apolys(j)
2006	polycoins(ip) = j
2007	ELSE IF ( dist == maxcoindist) THEN
2008	IF (apolys(j) > maxcoinarea) THEN
2009	polycoins(ip) = j
2010	maxcoinarea = apolys(j)
2011	END IF
2012	END IF
2013	END IF
2014	END DO
2015	END IF
2016	END DO
2017
2018	RETURN
2019
2020	END SUBROUTINE coincidence_all_polys
2021
2022	SUBROUTINE coincidence_poly(dbg, dx, dy, poly, Npoly, polys, polycoin, coinNpts)
2023	! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
2024
2025	IMPLICIT NONE
2026
2027	LOGICAL, INTENT(in) :: dbg
2028	INTEGER, INTENT(in) :: dx, dy, Npoly
2029	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: poly
2030	INTEGER, DIMENSION(dx,dy), INTENT(in) :: polys
2031	INTEGER, INTENT(out) :: polycoin
2032	INTEGER, DIMENSION(Npoly), INTENT(out) :: coinNpts
2033
2034	! Local
2035	INTEGER :: i, j, ip
2036	INTEGER :: maxcorr
2037	INTEGER :: Nmaxcorr
2038
2039	!!!!!!! Variables
2040	! dx,dy: dimension of the space
2041	! poly: bolean polygon to meet
2042	! Npoly: number of polygons on the 2D space
2043	! polys: 2D field of polygons identified by their integer number (0 for no polygon)
2044	! polycoin: coincident polyogn
2045	! -1: no-coincidence
2046	! 1 < Npoly: single coincidence with a given polygon
2047	! -9: coincidence with more than one polygon
2048	! coinNpts: number of points coincident with each polygon
2049
2050	fname = 'coincidence_poly'
2051	IF (dbg) PRINT *,TRIM(fname)
2052
2053	IF (dbg) THEN
2054	PRINT *,' Boolean polygon to search coincidences ...'
2055	DO i=1,dx
2056	PRINT *,poly(i,:)
2057	END DO
2058
2059	PRINT *,' 2D polygons space ...'
2060	DO i=1,dx
2061	PRINT '(1000(I7,1x))',polys(i,:)
2062	END DO
2063	END IF
2064
2065	! Looking for coincient points for the polygon
2066	coinNpts = 0
2067	DO i=1,dx
2068	DO j=1,dy
2069	IF (poly(i,j) .AND. polys(i,j) .NE. 0) coinNpts(polys(i,j)) = coinNpts(polys(i,j)) + 1
2070	END DO
2071	END DO
2072
2073	maxcorr = 0
2074	maxcorr = INT(MAXVAL(coinNpts*1.))
2075
2076	IF (dbg) PRINT *,' Maximum coincidence:', maxcorr
2077	IF (maxcorr == 0) THEN
2078	polycoin = -1
2079	ELSE
2080	Nmaxcorr = 0
2081	DO ip=1, Npoly
2082	IF (coinNpts(ip) == maxcorr) THEN
2083	Nmaxcorr=Nmaxcorr+1
2084	polycoin = ip
2085	END IF
2086	END DO
2087	IF (Nmaxcorr > 1) polycoin = -9
2088	END IF
2089
2090	IF (dbg) THEN
2091	PRINT *,' Coincidences for each polygon _______'
2092	DO ip=1, Npoly
2093	PRINT *,' ', ip,': ', coinNpts(ip)
2094	END DO
2095	END IF
2096
2097	RETURN
2098
2099	END SUBROUTINE coincidence_poly
2100
2101	SUBROUTINE all_polygons_properties(dbg, dx, dy, Npoly, polys, lon, lat, values, xres, yres, projN, &
2102	Npolyptss, xxtrms, yxtrms, meanctrs, meanwctrs, areas, nvals, xvals, mvals, m2vals, stdvals, &
2103	Nquant, quants, nvcoords, xvcoords, meanvnctrs, meanvxctrs)
2104	! Subroutine to determine the properties of all polygons in a 2D field:
2105	! Number of grid points
2106	! grid-point coordinates of the minimum and maximum of the path along x,y axes
2107	! grid coordinates of center from the mean of the coordinates of the poygon locations
2108	! lon, lat center from the area weighted mean of the coordinates of the polygon locations
2109	! area of the polygon (km2)
2110	! minimum and maximum of the values within the polygon
2111	! mean of the values within the polygon
2112	! quadratic mean of the values within the polygon
2113	! standard deviation of the values within the polygon
2114	! number of quantiles
2115	! quantiles of the values within the polygon
2116	! grid coordinates of the minimum, maximum value within the polygon
2117	! lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon
2118
2119	IMPLICIT NONE
2120
2121	LOGICAL, INTENT(in) :: dbg
2122	INTEGER, INTENT(in) :: dx, dy, Npoly, Nquant
2123	INTEGER, DIMENSION(dx,dy), INTENT(in) :: polys
2124	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: lon, lat, values
2125	REAL(r_k), INTENT(in) :: xres, yres
2126	CHARACTER(len=1000), INTENT(in) :: projN
2127	INTEGER, DIMENSION(Npoly), INTENT(out) :: Npolyptss
2128	INTEGER, DIMENSION(Npoly,2), INTENT(out) :: xxtrms, yxtrms, meanctrs
2129	REAL(r_k), DIMENSION(Npoly), INTENT(out) :: areas
2130	REAL(r_k), DIMENSION(Npoly), INTENT(out) :: nvals, xvals, mvals, m2vals, stdvals
2131	REAL(r_k), DIMENSION(Npoly, Nquant), INTENT(out) :: quants
2132	INTEGER, DIMENSION(Npoly,2), INTENT(out) :: nvcoords, xvcoords
2133	REAL(r_k), DIMENSION(Npoly,2), INTENT(out) :: meanwctrs, meanvnctrs, meanvxctrs
2134
2135	! Local
2136	INTEGER :: ip
2137	LOGICAL, DIMENSION(dx,dy) :: boolpoly
2138
2139	!!!!!!! Variables
2140	! dx,dy: size of the space
2141	! Npoly: number of polygons
2142	! polys: polygon matrix with all polygons (as integer number per polygon)
2143	! lon, lat: geographical coordinates of the grid-points of the matrix
2144	! values: values of the 2D field to use
2145	! [x/y]res resolution along the x and y axis
2146	! projN: name of the projection
2147	! 'ctsarea': Constant Area
2148	! 'lon/lat': for regular longitude-latitude
2149	! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
2150	! Npolyptss: number of points of the polygons
2151	! [x/y]xtrms: grid-point coordinates of minimum and maximum coordinates of the polygons
2152	! meanctrs: center from the mean of the coordinates of the polygons
2153	! meanwctrs: lon, lat coordinates of the center from the spatial-weighted mean of the polygons
2154	! areas: area of the polygons [km]
2155	! [n/x]vals: minimum and maximum of the values within the polygons
2156	! mvals: mean of the values within the polygons
2157	! m2vals: quadratic mean of the values within the polygons
2158	! stdvals: standard deviation of the values within the polygons
2159	! Nquant: number of quantiles
2160	! quants: quantiles of the values within the polygons
2161	! [n/x]vcoords: grid coordinates of the minimum/maximum of the values within the polygons
2162	! meanv[n/x]ctrs: lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygons
2163
2164	fname = 'all_polygons_properties'
2165
2166	! Initializing matrices
2167	Npolyptss = -1
2168	xxtrms = fillval64
2169	yxtrms = fillval64
2170	meanctrs = fillval64
2171	meanwctrs = fillval64
2172	areas = fillval64
2173	nvals = fillvalI
2174	xvals = fillval64
2175	mvals = fillval64
2176	m2vals = fillval64
2177	stdvals = fillval64
2178	quants = fillval64
2179	nvcoords = fillvalI
2180	xvcoords = fillvalI
2181	meanvnctrs = fillval64
2182	meanvxctrs = fillval64
2183
2184	DO ip=1, Npoly
2185	boolpoly = polys == ip
2186	CALL polygon_properties(dbg, dx, dy, boolpoly, lon, lat, values, xres, yres, projN, Npolyptss(ip),&
2187	xxtrms(ip,:), yxtrms(ip,:), meanctrs(ip,:), meanwctrs(ip,:), areas(ip), nvals(ip), xvals(ip), &
2188	mvals(ip), m2vals(ip), stdvals(ip), Nquant, quants(ip,:), nvcoords(ip,:), xvcoords(ip,:), &
2189	meanvnctrs(ip,:), meanvxctrs(ip,:))
2190	END DO
2191
2192	RETURN
2193
2194	END SUBROUTINE all_polygons_properties
2195
2196	SUBROUTINE polygon_properties(dbg, dx, dy, poly, lon, lat, values, xres, yres, projN, Npolypts, &
2197	xxtrm, yxtrm, meanctr, meanwctr, area, nval, xval, mval, m2val, stdval, Nquant, quant, nvcoord, &
2198	xvcoord, meanvnctr, meanvxctr)
2199	! Subroutine to determine the properties of a polygon (as .TRUE. matrix)
2200	! Number of grid points
2201	! grid-point coordinates of the minimum and maximum of the path along x,y axes
2202	! grid coordinates of center from the mean of the coordinates of the poygon locations
2203	! lon, lat center from the area weighted mean of the coordinates of the polygon locations
2204	! area of the polygon (km2)
2205	! minimum and maximum of the values within the polygon
2206	! mean of the values within the polygon
2207	! quadratic mean of the values within the polygon
2208	! standard deviation of the values within the polygon
2209	! number of quantiles
2210	! quantiles of the values within the polygon
2211	! grid coordinates of the minimum, maximum value within the polygon
2212	! lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon
2213
2214	IMPLICIT NONE
2215
2216	LOGICAL, INTENT(in) :: dbg
2217	INTEGER, INTENT(in) :: dx, dy, Nquant
2218	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: poly
2219	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: lon, lat, values
2220	REAL(r_k), INTENT(in) :: xres, yres
2221	CHARACTER(len=1000), INTENT(in) :: projN
2222	INTEGER, INTENT(out) :: Npolypts
2223	INTEGER, DIMENSION(2), INTENT(out) :: xxtrm, yxtrm, meanctr
2224	INTEGER, DIMENSION(2), INTENT(out) :: nvcoord, xvcoord
2225	REAL(r_k), DIMENSION(2), INTENT(out) :: meanwctr, meanvnctr, meanvxctr
2226	REAL(r_k), INTENT(out) :: area
2227	REAL(r_k), INTENT(out) :: nval, xval, mval, m2val, stdval
2228	REAL(r_k), DIMENSION(Nquant), INTENT(out) :: quant
2229
2230	! Local
2231	INTEGER :: i, j, ip
2232	INTEGER :: ierr
2233	INTEGER, DIMENSION(:,:), ALLOCATABLE :: polypts
2234	REAL(r_k), DIMENSION(:), ALLOCATABLE :: polyvals, distvn, distvx
2235	REAL(r_k) :: lonNADIR, latNADIR
2236	REAL(r_k) :: sumRESx, sumRESy
2237	REAL(r_k), DIMENSION(dx,dy) :: xcorr, ycorr
2238	CHARACTER(len=200), DIMENSION(3) :: satSvals
2239	CHARACTER(len=50) :: projS
2240	REAL(r_k) :: sumDISTnlon, sumDISTnlat, sumDISTxlon, &
2241	sumDISTxlat
2242
2243	!!!!!!! Variables
2244	! dx,dy: size of the space
2245	! poly: polygon matrix (boolean)
2246	! lon, lat: geographical coordinates of the grid-points of the matrix
2247	! values: values of the 2D field to use
2248	! [x/y]res resolution along the x and y axis
2249	! projN: name of the projection
2250	! 'ctsarea': Constant Area
2251	! 'lon/lat': for regular longitude-latitude
2252	! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
2253	! Npolypts: number of points of the polygon
2254	! [x/y]xtrm: grid-point coordinates of minimum and maximum coordinates of the polygon
2255	! meanctr: center from the mean of the coordinates of the polygon
2256	! meanwctr: lon, lat coordinates of the center from the spatial-weighted mean of the polygon
2257	! area: area of the polygon [km]
2258	! [n/x]val: minimum and maximum of the values within the polygon
2259	! mval: mean of the values within the polygon
2260	! m2val: quadratic mean of the values within the polygon
2261	! stdval: standard deviation of the values within the polygon
2262	! Nquant: number of quantiles
2263	! quant: quantiles of the values within the polygon
2264	! [n/x]vcoord: grid coordinates of the minimum/maximum of the values within the polygon
2265	! meanv[n/x]ctr: lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon
2266
2267	fname = 'polygon_properties'
2268
2269	IF (dbg) PRINT *," '" // TRIM(fname) // "' ..."
2270
2271	! Getting grid-point coordinates of the polygon
2272	Npolypts = COUNT(poly)
2273
2274	IF (ALLOCATED(polypts)) DEALLOCATE(polypts)
2275	ALLOCATE(polypts(Npolypts,2), STAT=ierr)
2276	msg = "Problems allocating 'polypts'"
2277	CALL ErrMsg(msg, fname, ierr)
2278
2279	IF (ALLOCATED(polyvals)) DEALLOCATE(polyvals)
2280	ALLOCATE(polyvals(Npolypts), STAT=ierr)
2281	msg = "Problems allocating 'polyvals'"
2282	CALL ErrMsg(msg, fname, ierr)
2283
2284	IF (ALLOCATED(distvn)) DEALLOCATE(distvn)
2285	ALLOCATE(distvn(Npolypts), STAT=ierr)
2286	msg = "Problems allocating 'distvn'"
2287	CALL ErrMsg(msg, fname, ierr)
2288
2289	IF (ALLOCATED(distvx)) DEALLOCATE(distvx)
2290	ALLOCATE(distvx(Npolypts), STAT=ierr)
2291	msg = "Problems allocating 'distvx'"
2292	CALL ErrMsg(msg, fname, ierr)
2293
2294	IF (projN(1:7) == 'lon/lat') THEN
2295	projS = projN
2296	ELSE IF (projN(1:7) == 'ctsarea') THEN
2297	projS = projN
2298	ELSE IF (projN(1:9) == 'nadir-sat') THEN
2299	projS = 'nadir-sat'
2300	CALL split(projN, ',', 3, satSvals)
2301	READ(satSvals(2),'(F200.100)')lonNadir
2302	READ(satSvals(3),'(F200.100)')latNadir
2303	IF (dbg) PRINT *," 'nadir-geostationary-satellite' based projection of data with nadir " // &
2304	"location at:", lonNadir, latNadir
2305	ELSE
2306	msg = "Projection '" // TRIM(projN) // "' not ready" // CHAR(10) // " available ones: " // &
2307	"'ctsarea', 'lon/lat', 'nadir-sat'"
2308	CALL ErrMsg(msg,fname,-1)
2309	END IF
2310
2311	area = 0.
2312	sumRESx = 0.
2313	sumRESy = 0.
2314	meanwctr = 0.
2315	meanvnctr = 0.
2316	meanvxctr = 0.
2317	xcorr = 0.
2318	ycorr = 0.
2319
2320	nval = fillval64
2321	xval = -fillval64
2322
2323	ip = 1
2324	DO i=1,dx
2325	DO j=1,dy
2326	IF (poly(i,j)) THEN
2327	polypts(ip,1) = i
2328	polypts(ip,2) = j
2329	polyvals(ip) = values(i,j)
2330	SELECT CASE (TRIM(projS))
2331	CASE ('ctsarea')
2332	! Constant Area
2333	xcorr(i,j) = 1.
2334	ycorr(i,j) = 1.
2335	CASE ('lon/lat')
2336	! Area as fixed yres and sinus-lat varying for xres
2337	! IF (KIND(xcorr(i,j)) == KIND(1.d0)) THEN
2338	! xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
2339	! ELSE
2340	xcorr(i,j) = ABS(SIN(lon(i,j)*DegRad))
2341	! END IF
2342	ycorr(i,j) = 1.
2343	CASE ('nadir-sat')
2344	! Area from nadir resolution and degrading as we get far from satellite's nadir
2345	! GOES-E: 0 N, 75 W
2346	! IF (KIND(xcorr(i,j)) == KIND(1.d0)) THEN
2347	! xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
2348	! ELSE
2349	xcorr(i,j) = ABS(SIN(lon(i,j)*DegRad))
2350	! END IF
2351	ycorr(i,j) = 1.
2352	END SELECT
2353	area = area + xresxcorr(i,j)yres*ycorr(i,j)
2354	meanwctr(1) = meanwctr(1) + lon(i,j)xresxcorr(i,j)
2355	meanwctr(2) = meanwctr(2) + lat(i,j)yresycorr(i,j)
2356	IF (nval > values(i,j)) THEN
2357	nvcoord(1) = i
2358	nvcoord(2) = j
2359	nval = values(i,j)
2360	END IF
2361	IF (xval < values(i,j)) THEN
2362	xvcoord(1) = i
2363	xvcoord(2) = j
2364	xval = values(i,j)
2365	END IF
2366	ip = ip + 1
2367	END IF
2368	END DO
2369	END DO
2370
2371	IF (dbg) THEN
2372	PRINT *,' grid_coord lon lat value _______ '
2373	DO ip=1, Npolypts
2374	PRINT *, polypts(ip,:), ';', lon(polypts(ip,1),polypts(ip,2)), lat(polypts(ip,1),polypts(ip,2)),&
2375	':', polyvals(ip)
2376	END DO
2377	END IF
2378
2379	sumRESx = xres*SUM(xcorr)
2380	sumRESy = yres*SUM(ycorr)
2381
2382	xxtrm = (/ MINVAL(polypts(:,1)), MAXVAL(polypts(:,1)) /)
2383	yxtrm = (/ MINVAL(polypts(:,2)), MAXVAL(polypts(:,2)) /)
2384	meanctr = (/ SUM(polypts(:,1))/Npolypts, SUM(polypts(:,2))/Npolypts /)
2385	meanwctr = (/ meanwctr(1)/sumRESx, meanwctr(2)/sumRESy /)
2386
2387	IF (dbg) THEN
2388	PRINT *,' mean grid center: ', meanctr, ' weighted mean center: ', meanwctr
2389	END IF
2390
2391	! Statistics of the values within the polygon
2392	CALL StatsR_K(Npolypts, polyvals, nval, xval, mval, m2val, stdval)
2393
2394	IF (dbg) THEN
2395	PRINT *,' minimum value: ', nval, ' maximum:', xval, ' mean:', mval
2396	PRINT *,' coor. minimum: ', nvcoord
2397	PRINT *,' coor. maximum: ', xvcoord
2398	END IF
2399
2400	! Mean center weighted to minimum and maximum values
2401	! IF (KIND(polyvals(1)) == KIND(1.d0)) THEN
2402	! distvn = DABS(polyvals - nval)
2403	! distvx = DABS(polyvals - xval)
2404	! ELSE
2405	distvn = ABS(polyvals - nval)
2406	distvx = ABS(polyvals - xval)
2407	! END IF
2408
2409	meanvnctr = 0.
2410	meanvxctr = 0.
2411	sumDISTnlon = 0.
2412	sumDISTnlat = 0.
2413	sumDISTxlon = 0.
2414	sumDISTxlat = 0.
2415
2416	ip = 1
2417	DO i=1,dx
2418	DO j=1,dy
2419	IF (poly(i,j)) THEN
2420	meanvnctr(1) = meanvnctr(1)+lon(i,j)distvn(ip)xres*xcorr(i,j)
2421	meanvnctr(2) = meanvnctr(2)+lat(i,j)distvn(ip)yres*ycorr(i,j)
2422
2423	meanvxctr(1) = meanvxctr(1)+lon(i,j)distvx(ip)xres*xcorr(i,j)
2424	meanvxctr(2) = meanvxctr(2)+lat(i,j)distvx(ip)yres*ycorr(i,j)
2425
2426	sumDISTnlon = sumDISTnlon + distvn(ip)xresxcorr(i,j)
2427	sumDISTnlat = sumDISTnlat + distvn(ip)yresycorr(i,j)
2428	sumDISTxlon = sumDISTxlon + distvx(ip)xresxcorr(i,j)
2429	sumDISTxlat = sumDISTxlat + distvx(ip)yresycorr(i,j)
2430
2431	ip = ip + 1
2432	END IF
2433	END DO
2434	END DO
2435
2436	meanvnctr = (/ meanvnctr(1)/sumDISTnlon, meanvnctr(2)/sumDISTnlat /)
2437	meanvxctr = (/ meanvxctr(1)/sumDISTxlon, meanvxctr(2)/sumDISTxlat /)
2438
2439	IF (dbg) THEN
2440	PRINT *,' mean center to minimum: ', meanvnctr, ' to maximum: ', meanvxctr
2441	END IF
2442
2443	! Quantiles of the values within the polygon
2444	quant = -9999.d0
2445	IF (Npolypts > Nquant) THEN
2446	CALL quantilesR_K(Npolypts, polyvals, Nquant, quant)
2447	ELSE
2448	CALL SortR_K(polyvals, Npolypts)
2449	END IF
2450
2451	DEALLOCATE (polypts)
2452	DEALLOCATE (polyvals)
2453
2454	RETURN
2455
2456	END SUBROUTINE polygon_properties
2457
2458	SUBROUTINE polygons_t(dbg, dx, dy, dt, boolmatt, polys, Npoly)
2459	! Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!
2460
2461	IMPLICIT NONE
2462
2463	INTEGER, INTENT(in) :: dx, dy, dt
2464	LOGICAL, DIMENSION(dx,dy,dt), INTENT(in) :: boolmatt
2465	LOGICAL, INTENT(in) :: dbg
2466	INTEGER, DIMENSION(dt), INTENT(out) :: Npoly
2467	INTEGER, DIMENSION(dx,dy,dt), INTENT(out) :: polys
2468
2469	! Local
2470	INTEGER :: i,it
2471
2472	!!!!!!! Variables
2473	! dx,dy: spatial dimensions of the space
2474	! boolmatt: boolean matrix tolook for the polygons (.TRUE. based)
2475	! polys: found polygons
2476	! Npoly: number of polygons found
2477
2478	fname = 'polygons'
2479
2480	IF (dbg) PRINT *,TRIM(fname)
2481
2482	polys = -1
2483	Npoly = 0
2484
2485	DO it=1,dt
2486	IF (ANY(boolmatt(:,:,it))) THEN
2487	IF (dbg) THEN
2488	PRINT *,' it:', it, ' num. TRUE:', COUNT(boolmatt(:,:,it)), 'bool _______'
2489	DO i=1,dx
2490	PRINT *,boolmatt(i,:,it)
2491	END DO
2492	END IF
2493	CALL polygons(dbg, dx, dy, boolmatt(:,:,it), polys(:,:,it), Npoly(it))
2494	ELSE
2495	IF (dbg) THEN
2496	PRINT *,' it:', it, " without '.TRUE.' values skipiing it!!"
2497	END IF
2498	END IF
2499	END DO
2500
2501	END SUBROUTINE polygons_t
2502
2503	SUBROUTINE polygons(dbg, dx, dy, boolmat, polys, Npoly)
2504	! Subroutine to search the polygons of a boolean field. FORTRAN based. 1st = 1!
2505
2506	IMPLICIT NONE
2507
2508	INTEGER, INTENT(in) :: dx, dy
2509	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: boolmat
2510	LOGICAL, INTENT(in) :: dbg
2511	INTEGER, INTENT(out) :: Npoly
2512	INTEGER, DIMENSION(dx,dy), INTENT(out) :: polys
2513
2514	! Local
2515	INTEGER :: i, j, ip, ipp, Nppt
2516	INTEGER :: ierr
2517	INTEGER, DIMENSION(:,:), ALLOCATABLE :: borders
2518	LOGICAL, DIMENSION(dx,dy) :: isborder, isbordery
2519	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: paths
2520	INTEGER :: Npath
2521	INTEGER, DIMENSION(:), ALLOCATABLE :: Nptpaths
2522	INTEGER, DIMENSION(2) :: xtrx, xtry, meanpth
2523	INTEGER :: Nvertx, Npts
2524	INTEGER, DIMENSION(:,:), ALLOCATABLE :: vertxs, points
2525	LOGICAL, DIMENSION(:), ALLOCATABLE :: isin
2526
2527	!!!!!!! Variables
2528	! dx,dy: spatial dimensions of the space
2529	! boolmat: boolean matrix tolook for the polygons (.TRUE. based)
2530	! polys: found polygons
2531	! Npoly: number of polygons found
2532
2533	fname = 'polygons'
2534
2535	polys = -1
2536
2537	! The mathematical maximum woiuld be dx*dy/4, but let's be optimistic... (sorry Jero)
2538	Nppt = dx*dy/10
2539
2540	IF (ALLOCATED(borders)) DEALLOCATE(borders)
2541	ALLOCATE(borders(Nppt,2), STAT=ierr)
2542	msg = "Problems allocating matrix 'borders'"
2543	CALL ErrMsg(msg, fname, ierr)
2544
2545	IF (ALLOCATED(paths)) DEALLOCATE(paths)
2546	ALLOCATE(paths(Nppt,Nppt,2), STAT=ierr)
2547	msg = "Problems allocating matrix 'paths'"
2548	CALL ErrMsg(msg, fname, ierr)
2549
2550	IF (ALLOCATED(Nptpaths)) DEALLOCATE(Nptpaths)
2551	ALLOCATE(Nptpaths(Nppt), STAT=ierr)
2552	msg = "Problems allocating matrix 'Nptpaths'"
2553	CALL ErrMsg(msg, fname, ierr)
2554
2555	! Filling with the points of all the space with .TRUE.
2556	Npts = COUNT(boolmat)
2557
2558	IF (ALLOCATED(points)) DEALLOCATE(points)
2559	ALLOCATE(points(Npts,2), STAT=ierr)
2560	msg = "Problems allocating matrix 'points'"
2561	CALL ErrMsg(msg, fname, ierr)
2562
2563	! We only want to localize that points 'inside'
2564	ip = 1
2565	DO i=1, dx
2566	DO j=1, dy
2567	IF (boolmat(i,j)) THEN
2568	points(ip,1) = i
2569	points(ip,2) = j
2570	ip = ip + 1
2571	END IF
2572	END DO
2573	END DO
2574
2575	CALL borders_matrixL(dx, dy, Nppt, boolmat, borders, isborder, isbordery)
2576	CALL paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)
2577
2578	Npoly = Npath
2579
2580	DO ip=1, Npath
2581	IF (ALLOCATED(vertxs)) DEALLOCATE(vertxs)
2582	ALLOCATE(vertxs(Nptpaths(ip),2))
2583	msg = "Problems allocating matrix 'vertxs'"
2584	CALL ErrMsg(msg, fname, ierr)
2585
2586	IF (ALLOCATED(isin)) DEALLOCATE(isin)
2587	ALLOCATE(isin(Npts), STAT=ierr)
2588	msg = "Problems allocating matrix 'isin'"
2589	CALL ErrMsg(msg, fname, ierr)
2590
2591	isin = .FALSE.
2592
2593	IF (dbg) PRINT *, ' path:', ip, ' N pts:', Nptpaths(ip)
2594
2595	CALL path_properties(dx, dy, boolmat, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), xtrx, xtry, &
2596	meanpth, 'y', Nvertx, vertxs)
2597
2598	IF (dbg) THEN
2599	PRINT *, ' properties _______'
2600	PRINT *, ' x-extremes:', xtrx
2601	PRINT *, ' y-extremes:', xtry
2602	PRINT *, ' center mean:', meanpth
2603	PRINT *, ' y-vertexs:', Nvertx,' ________'
2604	DO i=1, Nvertx
2605	PRINT *,' ',i,':',vertxs(i,:)
2606	END DO
2607	END IF
2608
2609	CALL gridpoints_InsidePolygon(dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), Nvertx,&
2610	xtrx, xtry, vertxs, Npts, points, isin)
2611
2612	! Filling polygons
2613	DO ipp=1, Npts
2614	IF (isin(ipp)) polys(points(ipp,1),points(ipp,2)) = ip
2615	END DO
2616
2617	IF (dbg) THEN
2618	PRINT *,' boolmat isborder isbordery polygon (',xtrx(1),',',xtry(1),')x(',xtrx(2),',',xtry(2), &
2619	') _______'
2620	DO i=xtrx(1), xtrx(2)
2621	PRINT *,i,':',boolmat(i,xtry(1):xtry(2)), ' border ', isborder(i,xtry(1):xtry(2)), &
2622	' isbordery ', isbordery(i,xtry(1):xtry(2)), ' polygon ', polys(i,xtry(1):xtry(2))
2623	END DO
2624	END IF
2625
2626	END DO
2627
2628	! Cleaning polygons matrix of not-used and paths around holes
2629	CALL clean_polygons(dx, dy, boolmat, polys, Npoly, dbg)
2630
2631	DEALLOCATE (borders)
2632	DEALLOCATE (Nptpaths)
2633	DEALLOCATE (paths)
2634	DEALLOCATE (vertxs)
2635	DEALLOCATE (points)
2636	DEALLOCATE (isin)
2637
2638	RETURN
2639
2640	END SUBROUTINE polygons
2641
2642	SUBROUTINE clean_polygons(dx, dy, Lmat, pols, Npols, dbg)
2643	! Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
2644
2645	IMPLICIT NONE
2646
2647	INTEGER, INTENT(in) :: dx, dy
2648	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat
2649	INTEGER, INTENT(inout) :: Npols
2650	INTEGER, DIMENSION(dx,dy), INTENT(inout) :: pols
2651	LOGICAL, INTENT(in) :: dbg
2652
2653	! Local
2654	INTEGER :: i,j,ip,iprm
2655	INTEGER, DIMENSION(Npols) :: origPol, NotPol, neigPol
2656	INTEGER :: ispol, NnotPol
2657	CHARACTER(len=4) :: ISa
2658
2659	!!!!!!! Variables
2660	! dx, dy: size of the space
2661	! Lmat: original bolean matrix from which the polygons come from
2662	! Npols: original number of polygons
2663	! pols: polygons space
2664
2665	fname = 'clean_polygons'
2666	IF (dbg) PRINT *," At '" // TRIM(fname) // "' ..."
2667
2668	origPol = -1
2669
2670	! Looking for polygons already in space
2671	NnotPol = 0
2672	DO ip=1, Npols
2673	ispol = COUNT(pols-ip == 0)
2674	IF (ispol > 0) THEN
2675	origPol(ip) = ip
2676	ELSE
2677	NnotPol = NnotPol + 1
2678	NotPol(NnotPol) = ip
2679	neigPol(NnotPol) = -1
2680	END IF
2681	END DO
2682
2683	IF (dbg) THEN
2684	PRINT *,' It should be:', Npols, ' polygons, but already there are:', Npols - NnotPol
2685	PRINT *,' Polygons to remove:', NotPol(1:NnotPol)
2686	END IF
2687
2688	! Looking for the hole border of a polygon. This is identify as such polygon point which along
2689	! y-axis has NpolygonA, Npolygon, .FALSE.
2690	DO i=1,dx
2691	DO j=2,dy-1
2692	IF ( (pols(i,j-1) /= pols(i,j) .AND. pols(i,j+1) == -1) .AND. (COUNT(NotPol-pols(i,j)==0)==0) &
2693	.AND. (pols(i,j) /= -1) .AND. (pols(i,j-1) /= -1)) THEN
2694	IF (dbg) PRINT *,' Polygon:', pols(i,j), ' to be removed at point (',i,',',j,'); j-1:', &
2695	pols(i,j-1), ' j:', pols(i,j), ' j+1:', pols(i,j+1)
2696	NnotPol = NnotPol + 1
2697	NotPol(NnotPol) = pols(i,j)
2698	neigPol(NnotPol) = pols(i,j-1)
2699	END IF
2700	END DO
2701	END DO
2702
2703	IF (dbg) THEN
2704	PRINT *,' It should be:', Npols, ' polygons, but already there are:', Npols - NnotPol
2705	PRINT *,' Polygons to remove after looking for fake border-of-hole polygons _______'
2706	DO i=1, NnotPol
2707	PRINT *, ' Polygon:', NotPol(i), ' to be replaced by:', neigPol(i)
2708	END DO
2709	END IF
2710
2711	! Removing polygons
2712	DO iprm=1, NnotPol
2713	IF (neigPol(iprm) == -1) THEN
2714	WHERE (pols == NotPol(iprm))
2715	pols = -1
2716	END WHERE
2717	IF (dbg) THEN
2718	PRINT *,' removing polygon:', NotPol(iprm)
2719	END IF
2720	ELSE
2721	WHERE (pols == NotPol(iprm))
2722	pols = neigPol(iprm)
2723	END WHERE
2724	IF (dbg) THEN
2725	PRINT *,' replacing polygon:', NotPol(iprm), ' by:', neigPol(iprm)
2726	END IF
2727	END IF
2728	END DO
2729
2730	! Re-numbering (descending values)
2731	DO i = 1, NnotPol
2732	iprm = MAXVAL(NotPol(1:NnotPol))
2733	WHERE(pols > iprm)
2734	pols = pols - 1
2735	END WHERE
2736	j = Index1DArrayI(NotPol, NnotPol, iprm)
2737	NotPol(j) = -9
2738	END DO
2739
2740	Npols = Npols - NnotPol
2741
2742	RETURN
2743
2744	END SUBROUTINE clean_polygons
2745
2746	SUBROUTINE path_properties(dx, dy, Lmat, Nptspth, pth, xxtrm, yxtrm, meanctr, axs, Nvrtx, vrtxs)
2747	! Subroutine to determine the properties of a path:
2748	! extremes: minimum and maximum of the path along x,y axes
2749	! meancenter: center from the mean of the coordinates of the paths locations
2750	! vertexs: path point, without neighbours along a given axis
2751
2752	IMPLICIT NONE
2753
2754	INTEGER, INTENT(in) :: dx, dy, Nptspth
2755	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat
2756	INTEGER, DIMENSION(Nptspth,2), INTENT(in) :: pth
2757	CHARACTER, INTENT(in) :: axs
2758	INTEGER, DIMENSION(2), INTENT(out) :: meanctr, xxtrm, yxtrm
2759	INTEGER, INTENT(out) :: Nvrtx
2760	INTEGER, DIMENSION(Nptspth,2), INTENT(out) :: vrtxs
2761
2762	! Local
2763	INTEGER :: i, ip, jp
2764	INTEGER :: neig1, neig2
2765
2766	!!!!!!! Variables
2767	! dx,dy: size of the space
2768	! Lmat: original matrix of logical values for the path
2769	! Nptspth: number of points of the path
2770	! pth: path coordinates (clockwise)
2771	! axs: axis of finding the vertex
2772	! [x/y]xtrm: minimum and maximum coordinates of the path
2773	! meanctr: center from the mean of the coordinates of the path
2774	! Nvrtx: Number of vertexs of the path
2775	! vrtxs: coordinates of the vertexs
2776
2777	fname = 'path_properties'
2778
2779	vrtxs = -1
2780	Nvrtx = 0
2781
2782	xxtrm = (/ MINVAL(pth(:,1)), MAXVAL(pth(:,1)) /)
2783	yxtrm = (/ MINVAL(pth(:,2)), MAXVAL(pth(:,2)) /)
2784	meanctr = (/ SUM(pth(:,1))/Nptspth, SUM(pth(:,2))/Nptspth /)
2785
2786	IF (axs == 'x' .OR. axs == 'X') THEN
2787	! Looking vertexs along x-axis
2788	DO i=1, Nptspth
2789	ip = pth(i,1)
2790	jp = pth(i,2)
2791	neig1 = 0
2792	neig2 = 0
2793	! W-point
2794	IF (ip == 1) THEN
2795	neig1 = -1
2796	ELSE
2797	IF (.NOT.Lmat(ip-1,jp)) neig1 = -1
2798	END IF
2799	! E-point
2800	IF (ip == dx) THEN
2801	neig2 = -1
2802	ELSE
2803	IF (.NOT.Lmat(ip+1,jp)) neig2 = -1
2804	END IF
2805
2806	IF (neig1 == -1 .AND. neig2 == -1) THEN
2807	Nvrtx = Nvrtx + 1
2808	vrtxs(Nvrtx,:) = (/ip,jp/)
2809	END IF
2810	END DO
2811	ELSE IF (axs == 'y' .OR. axs == 'Y') THEN
2812	! Looking vertexs along x-axis
2813	DO i=1, Nptspth
2814	ip = pth(i,1)
2815	jp = pth(i,2)
2816
2817	neig1 = 0
2818	neig2 = 0
2819	! S-point
2820	IF (jp == 1) THEN
2821	neig1 = -1
2822	ELSE
2823	IF (.NOT.Lmat(ip,jp-1)) neig1 = -1
2824	END IF
2825	! N-point
2826	IF (jp == dy) THEN
2827	neig2 = -1
2828	ELSE
2829	IF (.NOT.Lmat(ip,jp+1)) neig2 = -1
2830	END IF
2831
2832	IF (neig1 == -1 .AND. neig2 == -1) THEN
2833	Nvrtx = Nvrtx + 1
2834	vrtxs(Nvrtx,:) = (/ ip, jp /)
2835	END IF
2836	END DO
2837	ELSE
2838	msg = "Axis '" // axs // "' not available" // CHAR(10) // " Available ones: 'x', 'X', 'y, 'Y'"
2839	CALL ErrMsg(msg, fname, -1)
2840	END IF
2841
2842	RETURN
2843
2844	END SUBROUTINE path_properties
2845
2846	SUBROUTINE gridpoints_InsidePolygon(dx, dy, isbrdr, Npath, path, Nvrtx, xpathxtrm, ypathxtrm, &
2847	vrtxs, Npts, pts, inside)
2848	! Subroutine to determine if a series of grid points are inside a polygon following ray casting algorithm
2849	! FROM: https://en.wikipedia.org/wiki/Point_in_polygon
2850
2851	IMPLICIT NONE
2852
2853	INTEGER, INTENT(in) :: dx,dy,Npath,Nvrtx,Npts
2854	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isbrdr
2855	INTEGER, DIMENSION(Npath,2), INTENT(in) :: path
2856	INTEGER, DIMENSION(2), INTENT(in) :: xpathxtrm, ypathxtrm
2857	INTEGER, DIMENSION(Npath,2) :: vrtxs
2858	INTEGER, DIMENSION(Npts,2), INTENT(in) :: pts
2859	LOGICAL, DIMENSION(Npts), INTENT(out) :: inside
2860
2861	! Local
2862	INTEGER :: i,j,ip,ix,iy
2863	INTEGER :: Nintersecs, isvertex, ispath
2864	INTEGER :: ierr
2865	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: halo_brdr
2866	INTEGER :: Nbrbrdr
2867
2868	!!!!!!! Variables
2869	! dx,dy: space size
2870	! Npath: number of points of the path of the polygon
2871	! path: path of the polygon
2872	! isbrdr: boolean matrix of the space wqith .T. on polygon border
2873	! Nvrtx: number of vertexs of the path
2874	! [x/y]pathxtrm extremes of the path
2875	! vrtxs: vertexs of the path along y-axis
2876	! Npts: number of points
2877	! pts: points to look for
2878	! inside: vector wether point is inside or not (coincident to a border is inside)
2879
2880	fname = 'gridpoints_InsidePolygon'
2881
2882	! Creation of a 1-grid point larger matrix to deal with points reaching the limits
2883	IF (ALLOCATED(halo_brdr)) DEALLOCATE(halo_brdr)
2884	ALLOCATE(halo_brdr(dx+2,dy+2), STAT=ierr)
2885	msg = "Problems allocating matrix 'halo_brdr'"
2886	CALL ErrMsg(msg, fname, ierr)
2887	halo_brdr = .FALSE.
2888
2889	DO i=1,dx
2890	halo_brdr(i+1,2:dy+1) = isbrdr(i,:)
2891	END DO
2892
2893	inside = .FALSE.
2894
2895	DO ip=1,Npts
2896	Nintersecs = 0
2897	ix = pts(ip,1)
2898	iy = pts(ip,2)
2899	! Point might be outside path range...
2900	IF (ix >= xpathxtrm(1) .AND. ix <= xpathxtrm(2) .AND. iy >= ypathxtrm(1) .AND. &
2901	iy <= ypathxtrm(2)) THEN
2902
2903	! It is a border point?
2904	ispath = index_list_coordsI(Npath, path, (/ix,iy/))
2905	IF (isbrdr(ix,iy) .AND. (ispath /= -1)) THEN
2906	inside(ip) = .TRUE.
2907	CYCLE
2908	END IF
2909
2910	! Looking along y-axis
2911	! Accounting for consecutives borders
2912	Nbrbrdr = 0
2913	DO j=MAX(1,ypathxtrm(1)-1),iy-1
2914	! Only counting that borders that are not vertexs
2915	ispath = index_list_coordsI(Npath, path, (/ix,j/))
2916	isvertex = index_list_coordsI(Npath, vrtxs, (/ix,j/))
2917
2918	IF (halo_brdr(ix+1,j+1) .AND. (ispath /= -1) .AND. (isvertex == -1) ) Nintersecs = Nintersecs + 1
2919	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
2920	Nbrbrdr = Nbrbrdr + 1
2921	ELSE
2922	! Will remove that consecutive borders above 2
2923	IF (Nbrbrdr /= 0) THEN
2924	Nintersecs = Nintersecs - MAX(Nbrbrdr-1, 0)
2925	Nbrbrdr = 0
2926	END IF
2927	END IF
2928	END DO
2929	IF (MOD(Nintersecs,2) /= 0) inside(ip) = .TRUE.
2930	END IF
2931
2932	END DO
2933
2934	RETURN
2935
2936	END SUBROUTINE gridpoints_InsidePolygon
2937
2938	SUBROUTINE look_clockwise_borders(dx,dy,Nbrdrs,brdrs,gbrdr,isbrdr,ix,iy,dbg,xf,yf,iff)
2939	! Subroutine to look clock-wise for a next point within a collection of borders (limits of a region)
2940
2941	IMPLICIT NONE
2942
2943	INTEGER, INTENT(in) :: dx, dy, Nbrdrs, ix, iy
2944	INTEGER, DIMENSION(Nbrdrs,2), INTENT(in) :: brdrs
2945	LOGICAL, DIMENSION(Nbrdrs), INTENT(in) :: gbrdr
2946	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isbrdr
2947	LOGICAL, INTENT(in) :: dbg
2948	INTEGER, INTENT(out) :: xf, yf, iff
2949
2950	! Local
2951	INTEGER :: isch
2952	CHARACTER(len=2), DIMENSION(8) :: Lclock
2953	INTEGER, DIMENSION(8,2) :: spt
2954	INTEGER :: iif, jjf
2955
2956	!!!!!!! Variables
2957	! dx, dy: 2D shape ot the space
2958	! Nbrdrs: number of brdrs found in this 2D space
2959	! brdrs: list of coordinates of the borders
2960	! gbrdr: accounts for the use if the given border point
2961	! isbrdr: accounts for the matrix of the point is a border or not
2962	! ix,iy: coordinates of the point to start to find for
2963	! xf,yf: coordinates of the found point
2964	! iff: position of the border found within the list of borders
2965
2966	fname = 'look_clockwise_borders'
2967
2968	! Looking clock-wise assuming that one starts from the westernmost point
2969
2970	! Label of the search
2971	lclock = (/ 'W ', 'NW', 'N ', 'NE', 'E ', 'SE', 'S ', 'SW' /)
2972	! Transformation to apply
2973	!spt = (/ (/-1,0/), (/-1,1/), (/0,1/), (/1,1/), (/1,0/), (/1,-1/), (/0,-1/), (/-1,-1/) /)
2974	spt(:,1) = (/ -1, -1, 0, 1, 1, 1, 0, -1 /)
2975	spt(:,2) = (/ 0, 1, 1, 1, 0, -1, -1, -1 /)
2976
2977	xf = -1
2978	yf = -1
2979	DO isch=1, 8
2980	! clock-wise search
2981	IF (spt(isch,1) >= 0) THEN
2982	iif = MIN(dx,ix+spt(isch,1))
2983	ELSE
2984	iif = MAX(1,ix+spt(isch,1))
2985	END IF
2986	IF (spt(isch,2) >= 0) THEN
2987	jjf = MIN(dy,iy+spt(isch,2))
2988	ELSE
2989	jjf = MAX(1,iy+spt(isch,2))
2990	END IF
2991	iff = index_list_coordsI(Nbrdrs, brdrs,(/iif,jjf/))
2992	IF (iff > 0) THEN
2993	IF (dbg) PRINT *,' ' // lclock(isch) // '-point:', iif,jjf, ':', iff, 'is',isbrdr(iif,jjf), &
2994	'got',gbrdr(iff)
2995	IF (isbrdr(iif,jjf) .AND. .NOT.gbrdr(iff)) THEN
2996	xf = iif
2997	yf = jjf
2998	EXIT
2999	END IF
3000	END IF
3001	END DO
3002
3003	RETURN
3004
3005	END SUBROUTINE look_clockwise_borders
3006
3007	SUBROUTINE borders_matrixL(dx,dy,dxy,Lmat,brdrs,isbrdr,isbrdry)
3008	! Subroutine to provide the borders of a logical array (interested in .TRUE.)
3009
3010	IMPLICIT NONE
3011
3012	INTEGER, INTENT(in) :: dx,dy,dxy
3013	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat
3014	INTEGER, DIMENSION(dxy,2), INTENT(out) :: brdrs
3015	LOGICAL, DIMENSION(dx,dy), INTENT(out) :: isbrdr, isbrdry
3016
3017	! Local
3018	INTEGER :: i,j,ib
3019
3020	!!!!!!! Variables
3021	! dx,dy: size of the space
3022	! dxy: maximum number of border points
3023	! Lmat: Matrix to look for the borders
3024	! brdrs: list of coordinates of the borders
3025	! isbrdr: matrix with .T./.F. wether the given matrix point is a border or not
3026	! isbrdry: matrix with .T./.F. wether the given matrix point is a border or not only along y-axis
3027
3028	fname = 'borders_matrixL'
3029
3030	isbrdr = .FALSE.
3031	brdrs = -1
3032	ib = 1
3033
3034	! Starting with the borders. If a given point is TRUE it is a path-vertex
3035	! Along y-axis
3036	DO i=1, dx
3037	IF (Lmat(i,1) .AND. .NOT.isbrdr(i,1)) THEN
3038	brdrs(ib,1) = i
3039	brdrs(ib,2) = 1
3040	isbrdr(i,1) = .TRUE.
3041	ib=ib+1
3042	END IF
3043	IF (Lmat(i,dy) .AND. .NOT.isbrdr(i,dy)) THEN
3044	brdrs(ib,1) = i
3045	brdrs(ib,2) = dy
3046	isbrdr(i,dy) = .TRUE.
3047	ib=ib+1
3048	END IF
3049	END DO
3050	! Along x-axis
3051	DO j=1, dy
3052	IF (Lmat(1,j) .AND. .NOT.isbrdr(1,j)) THEN
3053	brdrs(ib,1) = 1
3054	brdrs(ib,2) = j
3055	isbrdr(1,j) = .TRUE.
3056	ib=ib+1
3057	END IF
3058	IF (Lmat(dx,j) .AND. .NOT.isbrdr(dx,j)) THEN
3059	brdrs(ib,1) = dx
3060	brdrs(ib,2) = j
3061	isbrdr(dx,j) = .TRUE.
3062	ib=ib+1
3063	END IF
3064	END DO
3065
3066	isbrdry = isbrdr
3067
3068	! Border as that when looking on x-axis points with Lmat(i) /= Lmat(i+1)
3069	DO i=1, dx-1
3070	DO j=1, dy-1
3071	IF ( Lmat(i,j) .NEQV. Lmat(i+1,j) ) THEN
3072	IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN
3073	brdrs(ib,1) = i
3074	brdrs(ib,2) = j
3075	isbrdr(i,j) = .TRUE.
3076	ib=ib+1
3077	ELSE IF (Lmat(i+1,j) .AND. .NOT.isbrdr(i+1,j)) THEN
3078	brdrs(ib,1) = i+1
3079	brdrs(ib,2) = j
3080	isbrdr(i+1,j) = .TRUE.
3081	ib=ib+1
3082	END IF
3083	END IF
3084	! y-axis
3085	IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN
3086	IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN
3087	brdrs(ib,1) = i
3088	brdrs(ib,2) = j
3089	isbrdr(i,j) = .TRUE.
3090	isbrdry(i,j) = .TRUE.
3091	ib=ib+1
3092	ELSE IF (Lmat(i,j+1) .AND. .NOT.isbrdr(i,j+1)) THEN
3093	brdrs(ib,1) = i
3094	brdrs(ib,2) = j+1
3095	isbrdr(i,j+1) = .TRUE.
3096	isbrdry(i,j+1) = .TRUE.
3097	ib=ib+1
3098	END IF
3099	END IF
3100	END DO
3101	END DO
3102
3103	DO i=1, dx-1
3104	DO j=1, dy-1
3105	! y-axis
3106	IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN
3107	IF (Lmat(i,j)) THEN
3108	isbrdry(i,j) = .TRUE.
3109	ELSE IF (Lmat(i,j+1)) THEN
3110	isbrdry(i,j+1) = .TRUE.
3111	END IF
3112	END IF
3113	END DO
3114	END DO
3115	! only y-axis adding bands of 2 grid points
3116	DO i=1, dx-1
3117	DO j=2, dy-2
3118	IF ( (Lmat(i,j) .EQV. Lmat(i,j+1)) .AND. (Lmat(i,j).NEQV.Lmat(i,j-1)) .AND. (Lmat(i,j).NEQV.Lmat(i,j+2)) ) THEN
3119	IF (Lmat(i,j)) THEN
3120	isbrdry(i,j) = .TRUE.
3121	isbrdry(i,j+1) = .TRUE.
3122	END IF
3123	END IF
3124	END DO
3125	END DO
3126
3127	RETURN
3128
3129	END SUBROUTINE borders_matrixL
3130
3131	SUBROUTINE paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)
3132	! Subroutine to search the paths of a border field.
3133
3134	IMPLICIT NONE
3135
3136	INTEGER, INTENT(in) :: dx, dy, Nppt
3137	LOGICAL, INTENT(in) :: dbg
3138	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isborder
3139	INTEGER, DIMENSION(Nppt,2), INTENT(in) :: borders
3140	INTEGER, DIMENSION(Nppt,Nppt,2), INTENT(out) :: paths
3141	INTEGER, INTENT(out) :: Npath
3142	INTEGER, DIMENSION(Nppt), INTENT(out) :: Nptpaths
3143
3144	! Local
3145	INTEGER :: i,j,k,ib
3146	INTEGER :: ierr
3147	INTEGER :: Nbrdr
3148	LOGICAL, DIMENSION(:), ALLOCATABLE :: gotbrdr, emptygotbrdr
3149	INTEGER :: iipth, ipath, ip, Nptspath
3150	INTEGER :: iib, jjb, iip, ijp, iif, jjf, iff
3151	LOGICAL :: found, finishedstep
3152
3153	!!!!!!! Variables
3154	! dx,dy: spatial dimensions of the space
3155	! Nppt: possible number of paths and points that the paths can have
3156	! isborder: boolean matrix which provide the borders of the polygon
3157	! borders: coordinates of the borders of the polygon
3158	! paths: coordinates of each found path
3159	! Npath: number of paths found
3160	! Nptpaths: number of points per path
3161
3162	fname = 'paths_border'
3163
3164	! Sarting matrix
3165	paths = -1
3166	Npath = 0
3167	Nptspath = 0
3168	Nptpaths = -1
3169
3170	ib=1
3171	finishedstep = .FALSE.
3172
3173	! Number of border points
3174	DO ib=1, Nppt
3175	IF (borders(ib,1) == -1 ) EXIT
3176	END DO
3177	Nbrdr = ib-1
3178
3179	IF (dbg) THEN
3180	PRINT *,' borders _______'
3181	DO i=1,Nbrdr
3182	PRINT *,' ',i,':',borders(i,:)
3183	END DO
3184	END IF
3185
3186	! Matrix which keeps track if a border point has been located
3187	IF (ALLOCATED(gotbrdr)) DEALLOCATE(gotbrdr)
3188	ALLOCATE(gotbrdr(Nbrdr), STAT=ierr)
3189	msg = "Problems allocating matrix 'gotbrdr'"
3190	CALL ErrMsg(msg, fname, ierr)
3191	IF (ALLOCATED(emptygotbrdr)) DEALLOCATE(emptygotbrdr)
3192	ALLOCATE(emptygotbrdr(Nbrdr), STAT=ierr)
3193	msg = "Problems allocating matrix 'emptygotbrdr'"
3194	CALL ErrMsg(msg, fname, ierr)
3195
3196	gotbrdr = .FALSE.
3197	emptygotbrdr = .FALSE.
3198
3199	! Starting the fun...
3200
3201	! Looking along the lines and when a border is found, starting from there in a clock-wise way
3202	iipth = 1
3203	ipath = 1
3204	DO ib=1,Nbrdr
3205	iib = borders(iipth,1)
3206	jjb = borders(iipth,2)
3207	! Starting new path
3208	newpath: IF (.NOT.gotbrdr(iipth)) THEN
3209	ip = 1
3210	Nptspath = 1
3211	paths(ipath,ip,:) = borders(iipth,:)
3212	gotbrdr(iipth) = .TRUE.
3213	! Looking for following clock-wise search
3214	! Not looking for W, because search starts from the W
3215	iip = iib
3216	ijp = jjb
3217	DO k=1,Nbrdr
3218	IF (dbg) PRINT *,ipath,'iip jip:', iip, ijp
3219	found = .FALSE.
3220	CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE.,iif,jjf,iff)
3221	IF (iif /= -1) THEN
3222	ip=ip+1
3223	paths(ipath,ip,:) = (/ iif,jjf /)
3224	found = .TRUE.
3225	gotbrdr(iff) = .TRUE.
3226	iip = iif
3227	ijp = jjf
3228	Nptspath = Nptspath + 1
3229	END IF
3230
3231	IF (dbg) THEN
3232	PRINT *,iib,jjb,' end of this round path:', ipath, '_____', gotbrdr
3233	DO i=1, Nptspath
3234	PRINT *,' ',i,':',paths(ipath,i,:)
3235	END DO
3236	END IF
3237	! If it is not found a next point, might be because it is a non-polygon related value
3238	IF (.NOT.found) THEN
3239	IF (dbg) PRINT *,'NOT FOUND !!!', gotbrdr
3240	! Are still there available borders?
3241	IF (ALL(gotbrdr) .EQV. .TRUE.) THEN
3242	finishedstep = .TRUE.
3243	Npath = ipath
3244	Nptpaths(ipath) = Nptspath
3245	EXIT
3246	ELSE
3247	Nptpaths(ipath) = Nptspath
3248	! Let's have a look if the previous points in the path have already some 'non-located' neighbourgs
3249	DO i=Nptspath,1,-1
3250	iip = paths(ipath,i,1)
3251	ijp = paths(ipath,i,2)
3252	CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE., iif, &
3253	jjf,iff)
3254	IF (iif /= -1 .AND. iff /= -1) THEN
3255	IF (dbg) PRINT *,' re-take path from point:', iif,',',jjf,' n-path:', iff
3256	found = .TRUE.
3257	iipth = index_list_coordsI(Nppt, borders, (/iip,ijp/))
3258	EXIT
3259	END IF
3260	END DO
3261	IF (.NOT.found) THEN
3262	! Looking for the next available border point for the new path
3263	DO i=1,Nbrdr
3264	IF (.NOT.gotbrdr(i)) THEN
3265	iipth = i
3266	EXIT
3267	END IF
3268	END DO
3269	IF (dbg) PRINT *,' Looking for next path starting at:', iipth, ' point:', &
3270	borders(iipth,:)
3271	ipath=ipath+1
3272	EXIT
3273	END IF
3274	END IF
3275	ELSE
3276	IF (dbg) PRINT *,' looking for next point...'
3277	END IF
3278	IF (finishedstep) EXIT
3279	END DO
3280	END IF newpath
3281	END DO
3282	Npath = ipath
3283	Nptpaths(ipath) = Nptspath
3284
3285	DEALLOCATE (gotbrdr)
3286	DEALLOCATE (emptygotbrdr)
3287
3288	RETURN
3289
3290	END SUBROUTINE paths_border
3291
3292	SUBROUTINE rand_sample(Nvals, Nsample, sample)
3293	! Subroutine to randomly sample a range of indices
3294
3295	IMPLICIT NONE
3296
3297	INTEGER, INTENT(in) :: Nvals, Nsample
3298	INTEGER, DIMENSION(Nsample), INTENT(out) :: sample
3299
3300	! Local
3301	INTEGER :: i, ind, jmax
3302	REAL, DIMENSION(Nsample) :: randv
3303	CHARACTER(len=50) :: fname
3304	LOGICAL :: found
3305	LOGICAL, DIMENSION(Nvals) :: issampled
3306	CHARACTER(len=256) :: msg
3307	CHARACTER(len=10) :: IS1, IS2
3308
3309	!!!!!!! Variables
3310	! Nvals: number of values
3311	! Nsamples: number of samples
3312	! sample: samnple
3313	fname = 'rand_sample'
3314
3315	IF (Nsample > Nvals) THEN
3316	WRITE(IS1,'(I10)')Nvals
3317	WRITE(IS2,'(I10)')Nsample
3318	msg = 'Sampling of ' // TRIM(IS1) // ' is too big for ' // TRIM(IS1) // 'values'
3319	CALL ErrMsg(msg, fname, -1)
3320	END IF
3321
3322	! Generation of random numbers always the same series during the whole program!
3323	CALL RANDOM_NUMBER(randv)
3324
3325	! Making sure that we do not repeat any value
3326	issampled = .FALSE.
3327
3328	DO i=1, Nsample
3329	! Generation of the index from the random numbers
3330	ind = MAX(INT(randv(i)*Nvals), 1)
3331
3332	IF (.NOT.issampled(ind)) THEN
3333	sample(i) = ind
3334	issampled(ind) = .TRUE.
3335	ELSE
3336	! Looking around the given index
3337	!PRINT *,' Index :', ind, ' already sampled!', issampled(ind)
3338	found = .FALSE.
3339	DO jmax=1, Nvals
3340	ind = MIN(ind+jmax, Nvals)
3341	IF (.NOT.issampled(ind)) THEN
3342	sample(i) = ind
3343	issampled(ind) = .TRUE.
3344	found = .TRUE.
3345	EXIT
3346	END IF
3347	ind = MAX(1, ind-jmax)
3348	IF (.NOT.issampled(ind)) THEN
3349	sample(i) = ind
3350	issampled(ind) = .TRUE.
3351	found = .TRUE.
3352	EXIT
3353	END IF
3354	END DO
3355	IF (.NOT.found) THEN
3356	msg = 'sampling could not be finished due to absence of available value!!'
3357	CALL ErrMsg(msg, fname, -1)
3358	END IF
3359	END IF
3360
3361	END DO
3362
3363	RETURN
3364
3365	END SUBROUTINE rand_sample
3366
3367	SUBROUTINE PrintQuantilesR_K(Nvals, vals, Nquants, qtvs, bspc)
3368	! Subroutine to print the quantiles of values REAL(r_k)
3369
3370	IMPLICIT NONE
3371
3372	INTEGER, INTENT(in) :: Nvals, Nquants
3373	REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals
3374	REAL(r_k), DIMENSION(Nquants), INTENT(in) :: qtvs
3375	CHARACTER(len=1000), OPTIONAL :: bspc
3376
3377	! Local
3378	INTEGER :: iq
3379	LOGICAL, DIMENSION(Nvals) :: search1, search2, search
3380	CHARACTER(len=6) :: RS1
3381	CHARACTER(len=50) :: fname
3382	CHARACTER(len=1000) :: bspcS
3383
3384	!!!!!!! Variables
3385	! vals: series of values
3386	! qtvs: values of the quantiles
3387	! bspc: base quantity of spaces
3388
3389	fname = 'PrintQuantilesR_K'
3390
3391	IF (PRESENT(bspc)) THEN
3392	bspcS = bspc
3393	ELSE
3394	bspcS = ' '
3395	END IF
3396
3397	DO iq=1, Nquants-1
3398
3399	WHERE (vals >= qtvs(iq))
3400	search1 = .TRUE.
3401	ELSEWHERE
3402	search1 = .FALSE.
3403	END WHERE
3404
3405	WHERE (vals < qtvs(iq+1))
3406	search2 = .TRUE.
3407	ELSEWHERE
3408	search2 = .FALSE.
3409	END WHERE
3410
3411	WHERE (search1 .AND. search2)
3412	search = .TRUE.
3413	ELSEWHERE
3414	search = .FALSE.
3415	END WHERE
3416
3417	WRITE(RS1, '(F6.2)')(iq)*100./(Nquants-1)
3418	PRINT *, TRIM(bspcS) // '[',iq,']', TRIM(RS1) // ' %:', qtvs(iq), 'N:', COUNT(search)
3419
3420	END DO
3421
3422	RETURN
3423
3424	END SUBROUTINE PrintQuantilesR_K
3425
3426	INTEGER FUNCTION FindMinimumR_K(x, dsize, Startv, Endv)
3427	! Function returns the location of the minimum in the section between Start and End.
3428
3429	IMPLICIT NONE
3430
3431	INTEGER, INTENT(in) :: dsize
3432	REAL(r_k), DIMENSION(dsize), INTENT(in) :: x
3433	INTEGER, INTENT(in) :: Startv, Endv
3434
3435	! Local
3436	REAL(r_k) :: Minimum
3437	INTEGER :: Location
3438	INTEGER :: i
3439
3440	Minimum = x(Startv) ! assume the first is the min
3441	Location = Startv ! record its position
3442	DO i = Startv+1, Endv ! start with next elements
3443	IF (x(i) < Minimum) THEN ! if x(i) less than the min?
3444	Minimum = x(i) ! Yes, a new minimum found
3445	Location = i ! record its position
3446	END IF
3447	END DO
3448
3449	FindMinimumR_K = Location ! return the position
3450
3451	END FUNCTION FindMinimumR_K
3452
3453	SUBROUTINE SwapR_K(a, b)
3454	! Subroutine swaps the values of its two formal arguments.
3455
3456	IMPLICIT NONE
3457
3458	REAL(r_k), INTENT(INOUT) :: a, b
3459	! Local
3460	REAL(r_k) :: Temp
3461
3462	Temp = a
3463	a = b
3464	b = Temp
3465
3466	END SUBROUTINE SwapR_K
3467
3468	SUBROUTINE SortR_K(x, Nx)
3469	! Subroutine receives an array x() r_K and sorts it into ascending order.
3470
3471	IMPLICIT NONE
3472
3473	INTEGER, INTENT(IN) :: Nx
3474	REAL(r_k), DIMENSION(Nx), INTENT(INOUT) :: x
3475
3476	! Local
3477	INTEGER :: i
3478	INTEGER :: Location
3479
3480	DO i = 1, Nx-1 ! except for the last
3481	Location = FindMinimumR_K(x, Nx-i+1, i, Nx) ! find min from this to last
3482	CALL SwapR_K(x(i), x(Location)) ! swap this and the minimum
3483	END DO
3484
3485	END SUBROUTINE SortR_K
3486
3487	SUBROUTINE quantilesR_K(Nvals, vals, Nquants, quants)
3488	! Subroutine to provide the quantiles of a given set of values of type real 'r_k'
3489
3490	IMPLICIT NONE
3491
3492	INTEGER, INTENT(in) :: Nvals, Nquants
3493	REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals
3494	REAL(r_k), DIMENSION(Nquants), INTENT(out) :: quants
3495
3496	! Local
3497	INTEGER :: i
3498	REAL(r_k) :: minv, maxv
3499	REAL(r_k), DIMENSION(Nvals) :: sortedvals
3500
3501	!!!!!!! Variables
3502	! Nvals: number of values
3503	! Rk: kind of real of the values
3504	! vals: values
3505	! Nquants: number of quants
3506	! quants: values at which the quantile start
3507
3508	minv = MINVAL(vals)
3509	maxv = MAXVAL(vals)
3510
3511	sortedvals = vals
3512	! Using from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
3513	CALL SortR_K(sortedvals, Nvals)
3514
3515	quants(1) = minv
3516	DO i=2, Nquants
3517	quants(i) = sortedvals(INT((i-1)*Nvals/Nquants))
3518	END DO
3519
3520	END SUBROUTINE quantilesR_K
3521
3522
3523	SUBROUTINE StatsR_K(Nvals, vals, minv, maxv, mean, mean2, stdev)
3524	! Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a
3525	! series of r_k numbers
3526
3527	IMPLICIT NONE
3528
3529	INTEGER, INTENT(in) :: Nvals
3530	REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals
3531	REAL(r_k), INTENT(out) :: minv, maxv, mean, mean2, stdev
3532
3533	!!!!!!! Variables
3534	! Nvals: number of values
3535	! vals: values
3536	! minv: minimum value of values
3537	! maxv: maximum value of values
3538	! mean: mean value of values
3539	! mean2: quadratic mean value of values
3540	! stdev: standard deviation of values
3541
3542	minv = MINVAL(vals)
3543	maxv = MAXVAL(vals)
3544
3545	mean=SUM(vals)
3546	mean2=SUM(vals*vals)
3547
3548	mean=mean/Nvals
3549	mean2=mean2/Nvals
3550
3551	stdev=SQRT(mean2 - mean*mean)
3552
3553	RETURN
3554
3555	END SUBROUTINE StatsR_k
3556
3557	END MODULE module_scientific

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format