Context Navigation

module_scientific.f90 @ 1784

Last change on this file since 1784 was 1782, checked in by lfita, 7 years ago
New version for matrix reconstruction taking as reference vector values
File size: 129.2 KB

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

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

Download in other formats:

Original Format