Context Navigation

module_scientific.f90 @ 1672

Last change on this file since 1672 was 1661, checked in by lfita, 8 years ago

Finally working version!

Three different methods have been added to determine the correct coordinate of the trajectory on multiple coincidencies 'mean': get coordinates from the areal-weighted mean of the centers of the given polygons and their areas 'largest': get the coorindates of the largest polygon 'closest': get the coordinates of the closest polygon

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

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

Download in other formats:

Original Format