Context Navigation

module_scientific.f90 @ 2074

Last change on this file since 2074 was 1913, checked in by lfita, 7 years ago

Adding:

`NcountR': Subroutine to count real values

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

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

Download in other formats:

Original Format