Context Navigation

module_scientific.f90 @ 1659

Last change on this file since 1659 was 1655, checked in by lfita, 8 years ago
Adding new structure with module_basic.f90 and infrastructure to deal with netCDF from FORTRAN
File size: 87.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	! quantilesR_K: Subroutine to provide the quantiles of a given set of values of type real 'r_k'
31	! rand_sample: Subroutine to randomly sample a range of indices
32	! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order.
33	! StatsR_K: Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a
34	! series of r_k numbers
35	! SwapR_K*: Subroutine swaps the values of its two formal arguments.
36
37	!!! *Functions/Subroutines to sort values adpated. The method used is usually referred to as "selection" method.
38	! from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
39
40	USE module_definitions
41	USE module_generic
42
43	CONTAINS
44
45	SUBROUTINE poly_overlap_tracks_area(dbg, dx, dy, dt, minarea, inNallpolys, allpolys, ctrpolys, &
46	areapolys, Nmaxpoly, Nmaxtracks, tracks, finaltracks)
47	! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum
48	! overlaping/coincidence filtrered by a minimal area
49
50	IMPLICIT NONE
51
52	LOGICAL, INTENT(in) :: dbg
53	INTEGER, INTENT(in) :: dx, dy, dt, Nmaxpoly, Nmaxtracks
54	INTEGER, DIMENSION(dt), INTENT(in) :: inNallpolys
55	INTEGER, DIMENSION(dx,dy,dt), INTENT(in) :: allpolys
56	REAL(r_k), INTENT(in) :: minarea
57	REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in) :: ctrpolys
58	REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in) :: areapolys
59	REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,dt), &
60	INTENT(out) :: tracks
61	REAL(r_k), DIMENSION(4,Nmaxtracks,dt), INTENT(out) :: finaltracks
62
63	! Local
64	INTEGER :: i, j, ip, it, iip, itt
65	INTEGER :: ierr
66	REAL(r_k), DIMENSION(Nmaxpoly) :: Aprevpolys, Acurrpolys
67	REAL(r_k), DIMENSION(2,Nmaxpoly) :: Cprevpolys, Ccurrpolys
68	INTEGER, DIMENSION(dt) :: Nallpolys
69	INTEGER, DIMENSION(dx,dy) :: meetpolys, searchpolys
70	INTEGER, DIMENSION(Nmaxpoly) :: coincidencies
71	INTEGER, DIMENSION(Nmaxpoly) :: prevID, currID
72	INTEGER, DIMENSION(:), ALLOCATABLE :: coins
73	INTEGER, DIMENSION(:,:), ALLOCATABLE :: coinsNpts
74	INTEGER :: Nmeet, Nsearch, Nindep
75	INTEGER, DIMENSION(dt) :: Nindeppolys
76	CHARACTER(len=5) :: NcoinS
77	INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt) :: polysIndep
78	INTEGER, DIMENSION(Nmaxpoly,dt) :: NpolysIndep
79	INTEGER, DIMENSION(Nmaxpoly,dt) :: SpolysIndep
80	INTEGER :: iindep, iiprev
81	INTEGER :: Nprev, NNprev, Ntprev
82	LOGICAL :: Indeppolychained
83	INTEGER :: itrack, ictrack
84	INTEGER :: ixp, iyp, ttrack
85	INTEGER, DIMENSION(dt) :: Ntracks
86	INTEGER :: idtrack, maxtrack
87
88	!!!!!!! Variables
89	! dx,dy,dt: space/time dimensions
90	! Nallpolys: Vector with the number of polygons at each time-step
91	! allpolys: Series of 2D field with the polygons
92	! minarea: minimal area (in same units as areapolys) to perform the tracking
93	! ctrpolys: center of the polygons
94	! areapolys: area of the polygons
95	! Nmaxpoly: Maximum possible number of polygons
96	! Nmaxtracks: maximum number of tracks
97	! tracks: series of consecutive polygons
98	! trackperiod: period of the track in time-steps
99
100	fname = 'poly_overlap_tracks_area'
101
102	IF (dbg) PRINT *,TRIM(fname)
103
104	! Number of independent polygons by time step
105	Nindeppolys = 0
106	! Number of polygons attached to each independent polygons by time step
107	NpolysIndep = 0
108	! ID of searching polygon attached to each independent polygons by time step
109	SpolysIndep = 0
110	! ID of polygons attached to each independent polygons by time step
111	polysIndep = 0
112	! ID of polygons from previous time-step
113	prevID = 0
114	! ID of polygons from current time-step
115	currID = 0
116
117	! First time-step all are independent polygons
118	it = 1
119	Nmeet = inNallpolys(it)
120	Nindeppolys(it) = Nmeet
121	ip = 0
122	meetpolys = allpolys(:,:,it)
123	DO i=1, Nmeet
124	IF (areapolys(i,it) >= minarea) THEN
125	ip = ip + 1
126	SpolysIndep(ip,it) = i
127	currID(ip) = i
128	Acurrpolys(ip) = areapolys(i,it)
129	Ccurrpolys(1,ip) = ctrpolys(1,i,it)
130	Ccurrpolys(2,ip) = ctrpolys(2,i,it)
131	NpolysIndep(ip,it) = 1
132	polysIndep(ip,1,it) = i
133	ELSE
134	WHERE (meetpolys == i)
135	meetpolys = 0
136	END WHERE
137	END IF
138	END DO
139	Nallpolys(1) = ip
140	Nindeppolys(1) = ip
141
142	! Starting step
143	it = 0
144	IF (dbg) THEN
145	PRINT *,' time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',0
146	PRINT *,' number of independent polygons:', Nindeppolys(it+1)
147	PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______'
148	DO i=1, Nindeppolys(it+1)
149	PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', &
150	polysIndep(i,1:NpolysIndep(i,it+1),it+1)
151	END DO
152	END IF
153
154	! Looking for the coincidencies at each time step
155	DO it=1, dt-1
156	! Number of times that a polygon has a coincidence
157	coincidencies = 0
158
159	Nmeet = inNallpolys(it+1)
160	searchpolys = meetpolys
161	meetpolys = allpolys(:,:,it+1)
162	prevID = 0
163	prevID = currID
164	Aprevpolys = Acurrpolys
165	Cprevpolys = Ccurrpolys
166	ip = 0
167
168	DO i=1, Nmeet
169	IF (areapolys(i,it+1) >= minarea) THEN
170	ip = ip + 1
171	currID(ip) = i
172	Acurrpolys(ip) = areapolys(i,it+1)
173	Acurrpolys(ip) = areapolys(i,it+1)
174	Ccurrpolys(1,ip) = ctrpolys(1,i,it+1)
175	Ccurrpolys(2,ip) = ctrpolys(2,i,it+1)
176	ELSE
177	WHERE (meetpolys == i)
178	meetpolys = 0
179	END WHERE
180	END IF
181	END DO
182	Nallpolys(it+1) = ip
183	Nindeppolys(it+1) = ip
184
185	Nmeet = Nallpolys(it+1)
186	! Looking throughout the independent polygons
187	Nsearch = Nindeppolys(it)
188
189	IF (ALLOCATED(coins)) DEALLOCATE(coins)
190	ALLOCATE(coins(Nmeet), STAT=ierr)
191	msg="Problems allocating 'coins'"
192	CALL ErrMsg(msg,fname,ierr)
193
194	IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts)
195	ALLOCATE(coinsNpts(Nmeet, Nsearch), STAT=ierr)
196	msg="Problems allocating 'coinsNpts'"
197	CALL ErrMsg(msg,fname,ierr)
198
199	CALL coincidence_all_polys_area(dbg, dx,dy, Nmeet, currID, meetpolys, Acurrpolys(1:Nmeet), &
200	Nsearch, prevID, searchpolys, Cprevpolys(:,1:Nsearch), Aprevpolys(1:Nsearch), coins, &
201	coinsNpts)
202
203	! Counting the number of times that a polygon has a coincidency
204	IF (dbg) THEN
205	PRINT *,' Coincidencies for the given time-step:', it+1,' _______'
206	DO i=1, Nmeet
207	PRINT *,currID(i), coins(i),' N search pts:', coinsNpts(i,:)
208	END DO
209	END IF
210
211	! Looking for the same equivalencies
212	Nindep = 0
213	DO i=1, Nmeet
214	IF (coins(i) == -1) THEN
215	Nindep = Nindep + 1
216	SpolysIndep(Nindep,it+1) = -1
217	NpolysIndep(Nindep,it+1) = NpolysIndep(Nindep,it+1) + 1
218	polysIndep(Nindep,NpolysIndep(Nindep,it+1),it+1) = currID(i)
219	ELSE IF (coins(i) == -9) THEN
220	WRITE(NcoinS,'(I5)')coins(i)
221	msg="coins= "//TRIM(NcoinS)//" This is an error. One should have always only one " // &
222	"coincidence of polygon"
223	CALL ErrMsg(msg, fname, -1)
224	ELSE
225	! Looking for coincidences with previous independent polygons
226	DO ip=1, Nsearch
227	! Looking into the polygons associated
228	NNprev = NpolysIndep(ip,it)
229	DO j=1, NNprev
230	IF (coins(i) == polysIndep(ip,j,it)) THEN
231	! Which index corresponds to this coincidence?
232	iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
233	IF (iindep == -1) THEN
234	Nindep = Nindep + 1
235	SpolysIndep(Nindep,it+1) = coins(i)
236	END IF
237	iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
238	IF (iindep < 0) THEN
239	PRINT *,' Looking for:', coins(i)
240	PRINT *,' Within:', SpolysIndep(1:Nindep,it+1)
241	PRINT *,' Might content:', polysIndep(ip,1:NNprev,it)
242	PRINT *,' From an initial list:', coins(1:Nmeet)
243	msg = 'Wrong index! There must be an index here'
244	CALL ErrMsg(msg,fname,iindep)
245	END IF
246	coincidencies(ip) = coincidencies(ip) + 1
247	NpolysIndep(iindep,it+1) = NpolysIndep(iindep,it+1) + 1
248	polysIndep(iindep,NpolysIndep(iindep,it+1),it+1) = currID(i)
249	EXIT
250	END IF
251	END DO
252	END DO
253	END IF
254	END DO
255	Nindeppolys(it+1) = Nindep
256
257	IF (dbg) THEN
258	PRINT *,' time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',Nsearch
259	PRINT *,' number of independent polygons:', Nindeppolys(it+1)
260	PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______'
261	DO i=1, Nindeppolys(it+1)
262	PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', &
263	polysIndep(i,1:NpolysIndep(i,it+1),it+1)
264	END DO
265	END IF
266	END DO
267
268	IF (dbg) THEN
269	PRINT *, 'Coincidencies to connect _______'
270	DO it=1, dt
271	PRINT *,' it:', it, ' Nindep:', Nindeppolys(it)
272	PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs'
273	DO ip=1, Nindeppolys(it)
274	PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':', &
275	polysIndep(ip,1:NpolysIndep(ip,it),it)
276	END DO
277	END DO
278
279	END IF
280
281	! Trajectories
282	! It should be done following the number of 'independent' polygons
283	! One would concatenate that independent polygons which share IDs from one step to another
284
285	! First time-step. Take all polygons
286	itrack = 0
287	tracks = 0.
288	Ntracks = 0
289	DO ip=1, Nindeppolys(1)
290	itrack = itrack + 1
291	tracks(1,1,itrack,1) = itrack*1.
292	tracks(2,1,itrack,1) = SpolysIndep(ip,1)
293	tracks(3,1,itrack,1) = ctrpolys(1,ip,1)
294	tracks(4,1,itrack,1) = ctrpolys(2,ip,1)
295	tracks(5,1,itrack,1) = 1
296	Ntracks(1) = Ntracks(1) + 1
297	END DO
298
299	! Looping allover already assigned tracks
300	timesteps: DO it=2, dt
301	IF (dbg) PRINT *,'track-timestep:', it, 'N indep polys:', Nindeppolys(it)
302	! Indep polygons current time-step
303	current_poly: DO i=1, Nindeppolys(it)
304	IF (dbg) PRINT *,' curent poly:', i, 'Prev poly:', SpolysIndep(i,it), ' N ass. polygons:', &
305	NpolysIndep(i,it), 'ass. poly:', polysIndep(i,1:NpolysIndep(i,it),it)
306	Indeppolychained = .FALSE.
307
308	! Number of tracks previous time-step
309	! Looping overall
310	it1_tracks: DO itt=1, Ntracks(it-1)
311	itrack = tracks(1,1,itt,it-1)
312	! Number polygons ID assigned
313	Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0)
314	IF (dbg) PRINT *,itt,' track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1)
315
316	! Looking for coincidencies
317	DO iip=1, Ntprev
318	IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN
319	Indeppolychained = .TRUE.
320	IF (dbg) PRINT *,' coincidence found by polygon:', tracks(2,iip,itt,it-1)
321	EXIT
322	END IF
323	END DO
324	IF (Indeppolychained) THEN
325	Ntracks(it) = Ntracks(it) + 1
326	ictrack = Ntracks(it)
327	! Assigning all the IDs to the next step of the track
328	DO iip=1, NpolysIndep(i,it)
329	iiprev = polysIndep(i,iip,it)
330	tracks(1,iip,ictrack,it) = itrack
331	tracks(2,iip,ictrack,it) = iiprev
332	ixp = ctrpolys(1,iiprev,it)
333	iyp = ctrpolys(2,iiprev,it)
334	tracks(3,iip,ictrack,it) = ixp
335	tracks(4,iip,ictrack,it) = iyp
336	tracks(5,iip,ictrack,it) = it
337	END DO
338	EXIT
339	END IF
340	IF (Indeppolychained) EXIT
341	END DO it1_tracks
342
343	! Creation of a new track
344	IF (.NOT.Indeppolychained) THEN
345	Ntracks(it) = Ntracks(it) + 1
346	ictrack = Ntracks(it)
347	! ID of new track
348	maxtrack = INT(MAXVAL(tracks(1,:,:,:)*1.))
349	IF (dbg) PRINT *,' New track!', maxtrack+1
350
351	! Assigning all the IDs to the next step of the track
352	DO j=1, NpolysIndep(i,it)
353	iiprev = polysIndep(i,j,it)
354	tracks(1,j,ictrack,it) = maxtrack+1
355	tracks(2,j,ictrack,it) = iiprev
356	ixp = ctrpolys(1,iiprev,it)
357	iyp = ctrpolys(2,iiprev,it)
358	tracks(3,j,ictrack,it) = ixp
359	tracks(4,j,ictrack,it) = iyp
360	tracks(5,j,ictrack,it) = it
361	END DO
362	END IF
363
364	END DO current_poly
365
366	IF (dbg) THEN
367	PRINT *,' At this time-step:', it, ' N trajectories:', Ntracks(it)
368	DO i=1, Ntracks(it)
369	Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0)
370	PRINT *,i ,'ID tracks:', tracks(1,1,i,it), 'ID polygon:', tracks(2,1:Nprev,i,it)
371	END DO
372	END IF
373
374	END DO timesteps
375
376	! Summarizing trajectories
377	! When multiple polygons are available, the mean of their central positions determines the position
378
379	finaltracks = 0.
380	maxtrack = MAXVAL(tracks(1,:,:,:))
381
382	DO it=1, dt
383	DO itt=1, Ntracks(it)
384	itrack = INT(tracks(1,1,itt,it))
385	Nprev = COUNT(INT(tracks(2,:,itt,it)) /= 0)
386	finaltracks(1,itrack,it) = itrack*1.
387	finaltracks(2,itrack,it) = SUM(tracks(3,:,itt,it))/Nprev*1.
388	finaltracks(3,itrack,it) = SUM(tracks(4,:,itt,it))/Nprev*1.
389	finaltracks(4,itrack,it) = it*1.
390	END DO
391	END DO
392
393	DEALLOCATE(coins)
394	DEALLOCATE(coinsNpts)
395
396	RETURN
397
398	END SUBROUTINE poly_overlap_tracks_area
399
400	SUBROUTINE coincidence_all_polys_area(dbg, dx, dy, Nallpoly, IDallpoly, allpoly, icpolys, Npoly, &
401	IDpolys, polys, cpolys, apolys, polycoins, coinNptss)
402	! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
403	! In case of multiple coincidencies, the closest and then the biggest is taken filtrered by a minimal area
404	! Here the difference is that the index does not coincide with its ID
405
406	IMPLICIT NONE
407
408	LOGICAL, INTENT(in) :: dbg
409	INTEGER, INTENT(in) :: dx, dy, Nallpoly, Npoly
410	INTEGER, DIMENSION(dx,dy), INTENT(in) :: allpoly, polys
411	INTEGER, DIMENSION(Nallpoly), INTENT(in) :: IDallpoly
412	INTEGER, DIMENSION(Npoly), INTENT(in) :: IDpolys
413	REAL(r_k), DIMENSION(2,Nallpoly), INTENT(in) :: icpolys
414	REAL(r_k), DIMENSION(2,Npoly), INTENT(in) :: cpolys
415	REAL(r_k), DIMENSION(Npoly), INTENT(in) :: apolys
416	INTEGER, DIMENSION(Nallpoly), INTENT(out) :: polycoins
417	INTEGER, DIMENSION(Nallpoly,Npoly), INTENT(out) :: coinNptss
418
419	! Local
420	INTEGER :: i, j, ip
421	INTEGER :: maxcorr
422	INTEGER :: Nmaxcorr
423	LOGICAL, DIMENSION(dx,dy) :: boolpoly
424	INTEGER :: maxcoin
425	REAL :: dist, maxcoindist, maxcoinarea
426	INTEGER, DIMENSION(Npoly) :: IDcoins
427
428	!!!!!!! Variables
429	! dx,dy: dimension of the space
430	! Nallpoly: Number of polygons to find coincidence
431	! allpoly: space with the polygons to meet
432	! icpolys: center of the polygons to look for the coincidence
433	! Npoly: number of polygons on the 2D space
434	! polys: 2D field of polygons identified by their integer number (0 for no polygon)
435	! cpolys: center of the polygons
436	! apolys: area of the polygons
437	! polycoins: coincident polyogn
438	! -1: no-coincidence
439	! 1 < Npoly: single coincidence with a given polygon
440	! -9: coincidence with more than one polygon
441	! coinNptss: number of points coincident with each polygon
442
443	fname = 'coincidence_all_polys_area'
444	IF (dbg) PRINT *,TRIM(fname)
445
446	DO ip=1, Nallpoly
447	boolpoly = allpoly == IDallpoly(ip)
448	CALL coincidence_poly_area(dbg, dx, dy, boolpoly, Npoly, polys, polycoins(ip), IDcoins, &
449	coinNptss(ip,:))
450	IF (dbg) PRINT *,' polygon', IDallpoly(ip), ' coincidence with:', polycoins(ip), 'IDpolys:', IDpolys(1:Npoly)
451
452	! Coincidence with more than one polygon
453	IF (polycoins(ip) == -9) THEN
454	maxcoindist = -10.
455	maxcoinarea = -10.
456	maxcoin = MAXVAL(coinNptss(ip,:))
457	DO j=1, Npoly
458	IF (coinNptss(ip,j) == maxcoin) THEN
459	dist = SQRT( (icpolys(1,ip)1.-cpolys(1,j)1.)*2 + (icpolys(2,ip)1.-cpolys(2,j)1.)*2 )
460	IF ( dist > maxcoindist) THEN
461	maxcoindist = dist
462	maxcoinarea = apolys(j)
463	polycoins(ip) = IDcoins(j)
464	ELSE IF ( dist == maxcoindist) THEN
465	IF (apolys(j) > maxcoinarea) THEN
466	polycoins(ip) = IDcoins(j)
467	maxcoinarea = apolys(j)
468	END IF
469	END IF
470	END IF
471	END DO
472	END IF
473	END DO
474
475	RETURN
476
477	END SUBROUTINE coincidence_all_polys_area
478
479	SUBROUTINE coincidence_poly_area(dbg, dx, dy, poly, Npoly, polys, polycoin, IDpoly, coinNpts)
480	! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
481	! Here the difference is that the index does not coincide with its ID
482
483	IMPLICIT NONE
484
485	LOGICAL, INTENT(in) :: dbg
486	INTEGER, INTENT(in) :: dx, dy, Npoly
487	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: poly
488	INTEGER, DIMENSION(dx,dy), INTENT(in) :: polys
489	INTEGER, INTENT(out) :: polycoin
490	INTEGER, DIMENSION(Npoly), INTENT(out) :: IDpoly, coinNpts
491
492	! Local
493	INTEGER :: i, j, ip
494	INTEGER :: maxcorr
495	INTEGER :: Nmaxcorr
496
497	!!!!!!! Variables
498	! dx,dy: dimension of the space
499	! poly: bolean polygon to meet
500	! Npoly: number of polygons on the 2D space
501	! polys: 2D field of polygons identified by their integer number (0 for no polygon)
502	! polycoin: coincident polyogn
503	! -1: no-coincidence
504	! 1 < Npoly: single coincidence with a given polygon
505	! -9: coincidence with more than one polygon
506	! IDpoly: ID of the found polygon
507	! coinNpts: number of points coincident with each polygon
508
509	fname = 'coincidence_poly_area'
510	IF (dbg) PRINT *,TRIM(fname)
511
512	IF (dbg) THEN
513	PRINT *,' Boolean polygon to search coincidences ...'
514	DO i=1,dx
515	PRINT *,poly(i,:)
516	END DO
517
518	PRINT *,' 2D polygons space ...'
519	DO i=1,dx
520	PRINT '(1000(I3,1x))',polys(i,:)
521	END DO
522	END IF
523
524	! Looking for coincient points for the polygon
525	coinNpts = 0
526	IDpoly = 0
527	ip = 0
528	DO i=1,dx
529	DO j=1,dy
530	IF (poly(i,j) .AND. polys(i,j) .NE. 0) THEN
531	IF (.NOT.ANY(IDpoly == polys(i,j))) THEN
532	ip = ip + 1
533	IDpoly(ip) = polys(i,j)
534	ELSE
535	ip = Index1DarrayI(IDpoly, Npoly, polys(i,j))
536	END IF
537	coinNpts(ip) = coinNpts(ip) + 1
538	END IF
539	END DO
540	END DO
541
542	maxcorr = 0
543	maxcorr = INT(MAXVAL(coinNpts*1.))
544
545	IF (dbg) PRINT *,' Maximum coincidence:', maxcorr
546	IF (maxcorr == 0) THEN
547	polycoin = -1
548	ELSE
549	Nmaxcorr = 0
550	DO ip=1, Npoly
551	IF (coinNpts(ip) == maxcorr) THEN
552	Nmaxcorr = Nmaxcorr+1
553	polycoin = IDpoly(ip)
554	END IF
555	END DO
556	IF (Nmaxcorr > 1) polycoin = -9
557	END IF
558
559	IF (dbg) THEN
560	PRINT *,' Coincidences for each polygon _______', Npoly
561	DO ip=1, Npoly
562	PRINT *,' ',ip, ' ID:', IDpoly(ip),': ', coinNpts(ip)
563	END DO
564	END IF
565
566	RETURN
567
568	END SUBROUTINE coincidence_poly_area
569
570	SUBROUTINE poly_overlap_tracks(dbg, dx, dy, dt, minarea, Nallpolys, allpolys, ctrpolys, &
571	areapolys, Nmaxpoly, Nmaxtracks, tracks, finaltracks)
572	! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum overlaping/coincidence
573
574	IMPLICIT NONE
575
576	LOGICAL, INTENT(in) :: dbg
577	INTEGER, INTENT(in) :: dx, dy, dt, Nmaxpoly, Nmaxtracks
578	INTEGER, DIMENSION(dt), INTENT(in) :: Nallpolys
579	INTEGER, DIMENSION(dx,dy,dt), INTENT(in) :: allpolys
580	REAL(r_k), INTENT(in) :: minarea
581	REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in) :: ctrpolys
582	REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in) :: areapolys
583	REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,dt), &
584	INTENT(out) :: tracks
585	REAL(r_k), DIMENSION(4,Nmaxtracks,dt), INTENT(out) :: finaltracks
586
587	! Local
588	INTEGER :: i, j, ip, it, iip, itt
589	INTEGER :: ierr
590	INTEGER, DIMENSION(Nmaxpoly,dt) :: coincidencies, NOcoincidencies
591	INTEGER, DIMENSION(:), ALLOCATABLE :: coins
592	INTEGER, DIMENSION(:,:), ALLOCATABLE :: coinsNpts
593	INTEGER, DIMENSION(Nmaxpoly,dt) :: polycoincidencies
594	INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt) :: coincidenciesNpts
595	INTEGER :: Nmeet, Nsearch, Nindep
596	INTEGER, DIMENSION(dt) :: Nindeppolys
597	CHARACTER(len=5) :: NcoinS
598	INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt) :: polysIndep
599	INTEGER, DIMENSION(Nmaxpoly,dt) :: NpolysIndep
600	INTEGER, DIMENSION(Nmaxpoly,dt) :: SpolysIndep
601	INTEGER :: iindep, iiprev
602	INTEGER :: Nprev, NNprev, Ntprev
603	LOGICAL :: Indeppolychained
604	INTEGER :: itrack, ictrack
605	INTEGER :: ixp, iyp, ttrack
606	INTEGER, DIMENSION(dt) :: Ntracks
607	INTEGER :: idtrack, maxtrack
608
609	!!!!!!! Variables
610	! dx,dy,dt: space/time dimensions
611	! Nallpolys: Vector with the number of polygons at each time-step
612	! allpolys: Series of 2D field with the polygons
613	! minarea: minimal area (in same units as areapolys) to perform the tracking
614	! ctrpolys: center of the polygons
615	! areapolys: area of the polygons
616	! Nmaxpoly: Maximum possible number of polygons
617	! Nmaxtracks: maximum number of tracks
618	! tracks: series of consecutive polygons
619	! trackperiod: period of the track in time-steps
620
621	fname = 'poly_overlap_tracks'
622
623	IF (dbg) PRINT *,TRIM(fname)
624
625	polycoincidencies = fillvalI
626	coincidenciesNpts = fillvalI
627	! Number of times that a polygon has a coincidence
628	coincidencies = 0
629	! Polygons without a coincidence
630	NOcoincidencies = 0
631	! Number of independent polygons by time step
632	Nindeppolys = 0
633	! Number of polygons attached to each independent polygons by time step
634	NpolysIndep = 0
635	! ID of searching polygon attached to each independent polygons by time step
636	SpolysIndep = 0
637	! ID of polygons attached to each independent polygons by time step
638	polysIndep = 0
639
640	! First time-step all are independent polygons
641	it = 1
642	Nmeet = Nallpolys(it)
643	Nindeppolys(it) = Nmeet
644	DO i=1, Nmeet
645	SpolysIndep(i,it) = i
646	NpolysIndep(1:Nmeet,it) = 1
647	polysIndep(1,i,it) = i
648	END DO
649
650	! Looking for the coincidencies at each time step
651	DO it=1, dt-1
652	Nmeet = Nallpolys(it+1)
653	Nsearch = Nallpolys(it)
654
655	IF (ALLOCATED(coins)) DEALLOCATE(coins)
656	ALLOCATE(coins(Nmeet), STAT=ierr)
657	msg="Problems allocating 'coins'"
658	CALL ErrMsg(msg,fname,ierr)
659
660	IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts)
661	ALLOCATE(coinsNpts(Nmeet, Nsearch), STAT=ierr)
662	msg="Problems allocating 'coinsNpts'"
663	CALL ErrMsg(msg,fname,ierr)
664
665	CALL coincidence_all_polys(dbg, dx, dy, Nmeet, allpolys(:,:,it+1), ctrpolys(:,1:Nmeet,it+1), &
666	Nsearch, allpolys(:,:,it), ctrpolys(:,1:Nsearch,it), areapolys(1:Nsearch,it), coins, coinsNpts)
667
668	polycoincidencies(1:Nmeet,it+1) = coins
669	coincidenciesNpts(1:Nmeet,1:Nsearch,it+1) = coinsNpts
670
671	! Counting the number of times that a polygon has a coincidency
672	IF (dbg) THEN
673	PRINT *,' Coincidencies for the given time-step:', it+1,' _______'
674	DO i=1, Nmeet
675	PRINT *,coins(i),' N search pts:', coinsNpts(i,:)
676	END DO
677	END IF
678
679	Nindep = 0
680	DO i=1, Nmeet
681	IF (coins(i) == -1) THEN
682	Nindep = Nindep + 1
683	NOcoincidencies(i,it+1) = 1
684	SpolysIndep(Nindep,it+1) = -1
685	NpolysIndep(Nindep,it+1) = NpolysIndep(Nindep,it+1) + 1
686	polysIndep(Nindep,NpolysIndep(Nindep,it+1),it+1) = i
687	ELSE IF (coins(i) == -9) THEN
688	WRITE(NcoinS,'(I5)')coins(i)
689	msg="coins= "//TRIM(NcoinS)//" This is an error. One should have always only one " // &
690	"coincidence of polygon"
691	CALL ErrMsg(msg, fname, -1)
692	ELSE
693	DO ip=1, Nsearch
694	IF (coins(i) == ip) THEN
695	IF (coincidencies(ip,it+1) == 0) THEN
696	Nindep = Nindep + 1
697	SpolysIndep(Nindep,it+1) = ip
698	END IF
699	coincidencies(ip,it+1) = coincidencies(ip,it+1) + 1
700	DO iindep=1, Nindep
701	IF (SpolysIndep(iindep,it+1) == coins(i)) THEN
702	NpolysIndep(iindep,it+1) = NpolysIndep(iindep,it+1) + 1
703	polysIndep(iindep,NpolysIndep(iindep,it+1),it+1) = i
704	END IF
705	END DO
706	END IF
707	END DO
708	END IF
709	END DO
710	Nindeppolys(it+1) = Nindep
711
712	IF (dbg) THEN
713	PRINT *,' time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',Nsearch
714	PRINT *,' number of independent polygons:', Nindeppolys(it+1)
715	PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______'
716	DO i=1, Nindeppolys(it+1)
717	PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', &
718	polysIndep(i,1:NpolysIndep(i,it+1),it+1)
719	END DO
720	END IF
721	END DO
722
723	IF (dbg) THEN
724	PRINT *, 'Coincidencies to connect _______'
725	DO it=1, dt
726	PRINT *,' it:', it, ' Nindep:', Nindeppolys(it)
727	PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs'
728	DO ip=1, Nindeppolys(it)
729	PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':', &
730	polysIndep(ip,1:NpolysIndep(ip,it),it)
731	END DO
732	END DO
733
734	END IF
735
736	! Trajectories
737	! It should be done following the number of 'independent' polygons
738	! One would concatenate that independent polygons which share IDs from one step to another
739
740	! First time-step. Take all polygons
741	itrack = 0
742	tracks = 0.
743	Ntracks = 0
744	DO ip=1, Nindeppolys(1)
745	itrack = itrack + 1
746	tracks(1,1,itrack,1) = itrack*1.
747	tracks(2,1,itrack,1) = SpolysIndep(ip,1)
748	tracks(3,1,itrack,1) = ctrpolys(1,ip,1)
749	tracks(4,1,itrack,1) = ctrpolys(2,ip,1)
750	tracks(5,1,itrack,1) = 1
751	Ntracks(1) = Ntracks(1) + 1
752	END DO
753
754	! Looping allover already assigned tracks
755	timesteps: DO it=2, dt
756	IF (dbg) PRINT *,'timestep:', it, 'N indep polys:', Nindeppolys(it)
757	! Indep polygons current time-step
758	current_poly: DO i=1, Nindeppolys(it)
759	IF (dbg) PRINT *,' curent poly:', i, 'Prev poly:', SpolysIndep(i,it), ' N ass. polygons:', &
760	NpolysIndep(i,it), 'ass. poly:', polysIndep(i,1:NpolysIndep(i,it),it)
761	Indeppolychained = .FALSE.
762
763	! Number of tracks previous time-step
764	! Looping overall
765	it1_tracks: DO itt=1, Ntracks(it-1)
766	itrack = tracks(1,1,itt,it-1)
767	! Number polygons ID assigned
768	Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0)
769	IF (dbg) PRINT *,itt,' track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1)
770
771	! Looking for coincidencies
772	DO iip=1, Ntprev
773	IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN
774	Indeppolychained = .TRUE.
775	IF (dbg) PRINT *,' coincidence found by polygon:', tracks(2,iip,itt,it-1)
776	EXIT
777	END IF
778	END DO
779	IF (Indeppolychained) THEN
780	Ntracks(it) = Ntracks(it) + 1
781	ictrack = Ntracks(it)
782	! Assigning all the IDs to the next step of the track
783	DO iip=1, NpolysIndep(i,it)
784	iiprev = polysIndep(i,iip,it)
785	tracks(1,iip,ictrack,it) = itrack
786	tracks(2,iip,ictrack,it) = iiprev
787	ixp = ctrpolys(1,iiprev,it)
788	iyp = ctrpolys(2,iiprev,it)
789	tracks(3,iip,ictrack,it) = ixp
790	tracks(4,iip,ictrack,it) = iyp
791	tracks(5,iip,ictrack,it) = it
792	END DO
793	EXIT
794	END IF
795	END DO it1_tracks
796
797	! Creation of a new track
798	IF (.NOT.Indeppolychained) THEN
799	Ntracks(it) = Ntracks(it) + 1
800	ictrack = Ntracks(it)
801	! ID of new track
802	maxtrack = INT(MAXVAL(tracks(1,:,:,:)*1.))
803	IF (dbg) PRINT *,' New track!', maxtrack+1
804
805	! Assigning all the IDs to the next step of the track
806	DO j=1, NpolysIndep(i,it)
807	iiprev = polysIndep(i,j,it)
808	tracks(1,j,ictrack,it) = maxtrack+1
809	tracks(2,j,ictrack,it) = iiprev
810	ixp = ctrpolys(1,iiprev,it)
811	iyp = ctrpolys(2,iiprev,it)
812	tracks(3,j,ictrack,it) = ixp
813	tracks(4,j,ictrack,it) = iyp
814	tracks(5,j,ictrack,it) = it
815	END DO
816	END IF
817
818	END DO current_poly
819
820	IF (dbg) THEN
821	PRINT *,' At this time-step:', it, ' N trajectories:', Ntracks(it)
822	DO i=1, Ntracks(it)
823	Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0)
824	PRINT *,i,'tracks:', tracks(2,1:Nprev,i,it)
825	END DO
826	END IF
827
828	END DO timesteps
829
830	! Summarizing trajectories
831	! When multiple polygons are available, the mean of their central positions determines the position
832
833	finaltracks = 0.
834	maxtrack = MAXVAL(tracks(1,:,:,:))
835
836	DO it=1, dt
837	DO itt=1, Ntracks(it)
838	itrack = INT(tracks(1,1,itt,it))
839	Nprev = COUNT(INT(tracks(2,:,itt,it)) /= 0)
840	PRINT *,'it:', it,'itrack:', itrack, 'Nprev:', Nprev
841	finaltracks(1,itrack,it) = itrack*1.
842	finaltracks(2,itrack,it) = SUM(tracks(3,:,itt,it))/Nprev
843	finaltracks(3,itrack,it) = SUM(tracks(4,:,itt,it))/Nprev
844	finaltracks(4,itrack,it) = it*1.
845	PRINT *,' finaltrack:', finaltracks(:,itrack,it)
846	END DO
847	END DO
848
849	RETURN
850
851	END SUBROUTINE poly_overlap_tracks
852
853	SUBROUTINE coincidence_all_polys(dbg, dx, dy, Nallpoly, allpoly, icpolys, Npoly, polys, cpolys, &
854	apolys, polycoins, coinNptss)
855	! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
856	! In case of multiple coincidencies, the closest and then the biggest is taken
857
858	IMPLICIT NONE
859
860	LOGICAL, INTENT(in) :: dbg
861	INTEGER, INTENT(in) :: dx, dy, Nallpoly, Npoly
862	INTEGER, DIMENSION(dx,dy), INTENT(in) :: allpoly, polys
863	REAL(r_k), DIMENSION(2,Nallpoly), INTENT(in) :: icpolys
864	REAL(r_k), DIMENSION(2,Npoly), INTENT(in) :: cpolys
865	REAL(r_k), DIMENSION(Npoly), INTENT(in) :: apolys
866	INTEGER, DIMENSION(Nallpoly), INTENT(out) :: polycoins
867	INTEGER, DIMENSION(Nallpoly,Npoly), INTENT(out) :: coinNptss
868
869	! Local
870	INTEGER :: i, j, ip
871	INTEGER :: maxcorr
872	INTEGER :: Nmaxcorr
873	LOGICAL, DIMENSION(dx,dy) :: boolpoly
874	INTEGER :: maxcoin
875	REAL :: dist, maxcoindist, maxcoinarea
876
877	!!!!!!! Variables
878	! dx,dy: dimension of the space
879	! Nallpoly: Number of polygons to find coincidence
880	! allpoly: space with the polygons to meet
881	! icpolys: center of the polygons to look for the coincidence
882	! Npoly: number of polygons on the 2D space
883	! polys: 2D field of polygons identified by their integer number (0 for no polygon)
884	! cpolys: center of the polygons
885	! apolys: area of the polygons
886	! polycoins: coincident polyogn
887	! -1: no-coincidence
888	! 1 < Npoly: single coincidence with a given polygon
889	! -9: coincidence with more than one polygon
890	! coinNptss: number of points coincident with each polygon
891
892	fname = 'coincidence_all_polys'
893	IF (dbg) PRINT *,TRIM(fname)
894
895	DO ip=1, Nallpoly
896	boolpoly = allpoly == ip
897	CALL coincidence_poly(dbg, dx, dy, boolpoly, Npoly, polys, polycoins(ip), coinNptss(ip,:))
898	IF (dbg) PRINT *,' polygon', ip, ' coincidence with:', polycoins(ip)
899
900	! Coincidence with more than one polygon
901	IF (polycoins(ip) == -9) THEN
902	maxcoindist = -10.
903	maxcoinarea = -10.
904	maxcoin = MAXVAL(coinNptss(ip,:))
905	DO j=1, Npoly
906	IF (coinNptss(ip,j) == maxcoin) THEN
907	dist = SQRT( (icpolys(1,ip)1.-cpolys(1,j)1.)*2 + (icpolys(2,ip)1.-cpolys(2,j)1.)*2 )
908	IF ( dist > maxcoindist) THEN
909	maxcoindist = dist
910	maxcoinarea = apolys(j)
911	polycoins(ip) = j
912	ELSE IF ( dist == maxcoindist) THEN
913	IF (apolys(j) > maxcoinarea) THEN
914	polycoins(ip) = j
915	maxcoinarea = apolys(j)
916	END IF
917	END IF
918	END IF
919	END DO
920	END IF
921	END DO
922
923	RETURN
924
925	END SUBROUTINE coincidence_all_polys
926
927	SUBROUTINE coincidence_poly(dbg, dx, dy, poly, Npoly, polys, polycoin, coinNpts)
928	! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons
929
930	IMPLICIT NONE
931
932	LOGICAL, INTENT(in) :: dbg
933	INTEGER, INTENT(in) :: dx, dy, Npoly
934	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: poly
935	INTEGER, DIMENSION(dx,dy), INTENT(in) :: polys
936	INTEGER, INTENT(out) :: polycoin
937	INTEGER, DIMENSION(Npoly), INTENT(out) :: coinNpts
938
939	! Local
940	INTEGER :: i, j, ip
941	INTEGER :: maxcorr
942	INTEGER :: Nmaxcorr
943
944	!!!!!!! Variables
945	! dx,dy: dimension of the space
946	! poly: bolean polygon to meet
947	! Npoly: number of polygons on the 2D space
948	! polys: 2D field of polygons identified by their integer number (0 for no polygon)
949	! polycoin: coincident polyogn
950	! -1: no-coincidence
951	! 1 < Npoly: single coincidence with a given polygon
952	! -9: coincidence with more than one polygon
953	! coinNpts: number of points coincident with each polygon
954
955	fname = 'coincidence_poly'
956	IF (dbg) PRINT *,TRIM(fname)
957
958	IF (dbg) THEN
959	PRINT *,' Boolean polygon to search coincidences ...'
960	DO i=1,dx
961	PRINT *,poly(i,:)
962	END DO
963
964	PRINT *,' 2D polygons space ...'
965	DO i=1,dx
966	PRINT '(1000(I3,1x))',polys(i,:)
967	END DO
968	END IF
969
970	! Looking for coincient points for the polygon
971	coinNpts = 0
972	DO i=1,dx
973	DO j=1,dy
974	IF (poly(i,j) .AND. polys(i,j) .NE. 0) coinNpts(polys(i,j)) = coinNpts(polys(i,j)) + 1
975	END DO
976	END DO
977
978	maxcorr = 0
979	maxcorr = INT(MAXVAL(coinNpts*1.))
980
981	IF (dbg) PRINT *,' Maximum coincidence:', maxcorr
982	IF (maxcorr == 0) THEN
983	polycoin = -1
984	ELSE
985	Nmaxcorr = 0
986	DO ip=1, Npoly
987	IF (coinNpts(ip) == maxcorr) THEN
988	Nmaxcorr=Nmaxcorr+1
989	polycoin = ip
990	END IF
991	END DO
992	IF (Nmaxcorr > 1) polycoin = -9
993	END IF
994
995	IF (dbg) THEN
996	PRINT *,' Coincidences for each polygon _______'
997	DO ip=1, Npoly
998	PRINT *,' ', ip,': ', coinNpts(ip)
999	END DO
1000	END IF
1001
1002	RETURN
1003
1004	END SUBROUTINE coincidence_poly
1005
1006	SUBROUTINE all_polygons_properties(dbg, dx, dy, Npoly, polys, lon, lat, values, xres, yres, projN, &
1007	Npolyptss, xxtrms, yxtrms, meanctrs, meanwctrs, areas, nvals, xvals, mvals, m2vals, stdvals, &
1008	Nquant, quants, nvcoords, xvcoords, meanvnctrs, meanvxctrs)
1009	! Subroutine to determine the properties of all polygons in a 2D field:
1010	! Number of grid points
1011	! grid-point coordinates of the minimum and maximum of the path along x,y axes
1012	! grid coordinates of center from the mean of the coordinates of the poygon locations
1013	! lon, lat center from the area weighted mean of the coordinates of the polygon locations
1014	! area of the polygon (km2)
1015	! minimum and maximum of the values within the polygon
1016	! mean of the values within the polygon
1017	! quadratic mean of the values within the polygon
1018	! standard deviation of the values within the polygon
1019	! number of quantiles
1020	! quantiles of the values within the polygon
1021	! grid coordinates of the minimum, maximum value within the polygon
1022	! lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon
1023
1024	IMPLICIT NONE
1025
1026	LOGICAL, INTENT(in) :: dbg
1027	INTEGER, INTENT(in) :: dx, dy, Npoly, Nquant
1028	INTEGER, DIMENSION(dx,dy), INTENT(in) :: polys
1029	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: lon, lat, values
1030	REAL(r_k), INTENT(in) :: xres, yres
1031	CHARACTER(len=1000), INTENT(in) :: projN
1032	INTEGER, DIMENSION(Npoly), INTENT(out) :: Npolyptss
1033	INTEGER, DIMENSION(Npoly,2), INTENT(out) :: xxtrms, yxtrms, meanctrs
1034	REAL(r_k), DIMENSION(Npoly), INTENT(out) :: areas
1035	REAL(r_k), DIMENSION(Npoly), INTENT(out) :: nvals, xvals, mvals, m2vals, stdvals
1036	REAL(r_k), DIMENSION(Npoly, Nquant), INTENT(out) :: quants
1037	INTEGER, DIMENSION(Npoly,2), INTENT(out) :: nvcoords, xvcoords
1038	REAL(r_k), DIMENSION(Npoly,2), INTENT(out) :: meanwctrs, meanvnctrs, meanvxctrs
1039
1040	! Local
1041	INTEGER :: ip
1042	LOGICAL, DIMENSION(dx,dy) :: boolpoly
1043
1044	!!!!!!! Variables
1045	! dx,dy: size of the space
1046	! Npoly: number of polygons
1047	! polys: polygon matrix with all polygons (as integer number per polygon)
1048	! lon, lat: geographical coordinates of the grid-points of the matrix
1049	! values: values of the 2D field to use
1050	! [x/y]res resolution along the x and y axis
1051	! projN: name of the projection
1052	! 'ctsarea': Constant Area
1053	! 'lon/lat': for regular longitude-latitude
1054	! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
1055	! Npolyptss: number of points of the polygons
1056	! [x/y]xtrms: grid-point coordinates of minimum and maximum coordinates of the polygons
1057	! meanctrs: center from the mean of the coordinates of the polygons
1058	! meanwctrs: lon, lat coordinates of the center from the spatial-weighted mean of the polygons
1059	! areas: area of the polygons [km]
1060	! [n/x]vals: minimum and maximum of the values within the polygons
1061	! mvals: mean of the values within the polygons
1062	! m2vals: quadratic mean of the values within the polygons
1063	! stdvals: standard deviation of the values within the polygons
1064	! Nquant: number of quantiles
1065	! quants: quantiles of the values within the polygons
1066	! [n/x]vcoords: grid coordinates of the minimum/maximum of the values within the polygons
1067	! 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
1068
1069	fname = 'all_polygons_properties'
1070
1071	! Initializing matrices
1072	Npolyptss = -1
1073	xxtrms = fillval64
1074	yxtrms = fillval64
1075	meanctrs = fillval64
1076	meanwctrs = fillval64
1077	areas = fillval64
1078	nvals = fillvalI
1079	xvals = fillval64
1080	mvals = fillval64
1081	m2vals = fillval64
1082	stdvals = fillval64
1083	quants = fillval64
1084	nvcoords = fillvalI
1085	xvcoords = fillvalI
1086	meanvnctrs = fillval64
1087	meanvxctrs = fillval64
1088
1089	DO ip=1, Npoly
1090	boolpoly = polys == ip
1091	CALL polygon_properties(dbg, dx, dy, boolpoly, lon, lat, values, xres, yres, projN, Npolyptss(ip),&
1092	xxtrms(ip,:), yxtrms(ip,:), meanctrs(ip,:), meanwctrs(ip,:), areas(ip), nvals(ip), xvals(ip), &
1093	mvals(ip), m2vals(ip), stdvals(ip), Nquant, quants(ip,:), nvcoords(ip,:), xvcoords(ip,:), &
1094	meanvnctrs(ip,:), meanvxctrs(ip,:))
1095	END DO
1096
1097	RETURN
1098
1099	END SUBROUTINE all_polygons_properties
1100
1101	SUBROUTINE polygon_properties(dbg, dx, dy, poly, lon, lat, values, xres, yres, projN, Npolypts, &
1102	xxtrm, yxtrm, meanctr, meanwctr, area, nval, xval, mval, m2val, stdval, Nquant, quant, nvcoord, &
1103	xvcoord, meanvnctr, meanvxctr)
1104	! Subroutine to determine the properties of a polygon (as .TRUE. matrix)
1105	! Number of grid points
1106	! grid-point coordinates of the minimum and maximum of the path along x,y axes
1107	! grid coordinates of center from the mean of the coordinates of the poygon locations
1108	! lon, lat center from the area weighted mean of the coordinates of the polygon locations
1109	! area of the polygon (km2)
1110	! minimum and maximum of the values within the polygon
1111	! mean of the values within the polygon
1112	! quadratic mean of the values within the polygon
1113	! standard deviation of the values within the polygon
1114	! number of quantiles
1115	! quantiles of the values within the polygon
1116	! grid coordinates of the minimum, maximum value within the polygon
1117	! lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon
1118
1119	IMPLICIT NONE
1120
1121	LOGICAL, INTENT(in) :: dbg
1122	INTEGER, INTENT(in) :: dx, dy, Nquant
1123	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: poly
1124	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: lon, lat, values
1125	REAL(r_k), INTENT(in) :: xres, yres
1126	CHARACTER(len=1000), INTENT(in) :: projN
1127	INTEGER, INTENT(out) :: Npolypts
1128	INTEGER, DIMENSION(2), INTENT(out) :: xxtrm, yxtrm, meanctr
1129	INTEGER, DIMENSION(2), INTENT(out) :: nvcoord, xvcoord
1130	REAL(r_k), DIMENSION(2), INTENT(out) :: meanwctr, meanvnctr, meanvxctr
1131	REAL(r_k), INTENT(out) :: area
1132	REAL(r_k), INTENT(out) :: nval, xval, mval, m2val, stdval
1133	REAL(r_k), DIMENSION(Nquant), INTENT(out) :: quant
1134
1135	! Local
1136	INTEGER :: i, j, ip
1137	INTEGER :: ierr
1138	INTEGER, DIMENSION(:,:), ALLOCATABLE :: polypts
1139	REAL(r_k), DIMENSION(:), ALLOCATABLE :: polyvals, distvn, distvx
1140	REAL(r_k) :: lonNADIR, latNADIR
1141	REAL(r_k) :: sumRESx, sumRESy
1142	REAL(r_k), DIMENSION(dx,dy) :: xcorr, ycorr
1143	CHARACTER(len=200), DIMENSION(3) :: satSvals
1144	CHARACTER(len=50) :: projS
1145	REAL(r_k) :: sumDISTnlon, sumDISTnlat, sumDISTxlon, &
1146	sumDISTxlat
1147
1148	!!!!!!! Variables
1149	! dx,dy: size of the space
1150	! poly: polygon matrix (boolean)
1151	! lon, lat: geographical coordinates of the grid-points of the matrix
1152	! values: values of the 2D field to use
1153	! [x/y]res resolution along the x and y axis
1154	! projN: name of the projection
1155	! 'ctsarea': Constant Area
1156	! 'lon/lat': for regular longitude-latitude
1157	! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
1158	! Npolypts: number of points of the polygon
1159	! [x/y]xtrm: grid-point coordinates of minimum and maximum coordinates of the polygon
1160	! meanctr: center from the mean of the coordinates of the polygon
1161	! meanwctr: lon, lat coordinates of the center from the spatial-weighted mean of the polygon
1162	! area: area of the polygon [km]
1163	! [n/x]val: minimum and maximum of the values within the polygon
1164	! mval: mean of the values within the polygon
1165	! m2val: quadratic mean of the values within the polygon
1166	! stdval: standard deviation of the values within the polygon
1167	! Nquant: number of quantiles
1168	! quant: quantiles of the values within the polygon
1169	! [n/x]vcoord: grid coordinates of the minimum/maximum of the values within the polygon
1170	! 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
1171
1172	fname = 'polygon_properties'
1173
1174	IF (dbg) PRINT *," '" // TRIM(fname) // "' ..."
1175
1176	! Getting grid-point coordinates of the polygon
1177	Npolypts = COUNT(poly)
1178
1179	IF (ALLOCATED(polypts)) DEALLOCATE(polypts)
1180	ALLOCATE(polypts(Npolypts,2), STAT=ierr)
1181	msg = "Problems allocating 'polypts'"
1182	CALL ErrMsg(msg, fname, ierr)
1183
1184	IF (ALLOCATED(polyvals)) DEALLOCATE(polyvals)
1185	ALLOCATE(polyvals(Npolypts), STAT=ierr)
1186	msg = "Problems allocating 'polyvals'"
1187	CALL ErrMsg(msg, fname, ierr)
1188
1189	IF (ALLOCATED(distvn)) DEALLOCATE(distvn)
1190	ALLOCATE(distvn(Npolypts), STAT=ierr)
1191	msg = "Problems allocating 'distvn'"
1192	CALL ErrMsg(msg, fname, ierr)
1193
1194	IF (ALLOCATED(distvx)) DEALLOCATE(distvx)
1195	ALLOCATE(distvx(Npolypts), STAT=ierr)
1196	msg = "Problems allocating 'distvx'"
1197	CALL ErrMsg(msg, fname, ierr)
1198
1199	IF (projN(1:7) == 'lon/lat') THEN
1200	projS = projN
1201	ELSE IF (projN(1:7) == 'ctsarea') THEN
1202	projS = projN
1203	ELSE IF (projN(1:9) == 'nadir-sat') THEN
1204	projS = 'nadir-sat'
1205	CALL split(projN, ',', 3, satSvals)
1206	READ(satSvals(2),'(F200.100)')lonNadir
1207	READ(satSvals(3),'(F200.100)')latNadir
1208	IF (dbg) PRINT *," 'nadir-geostationary-satellite' based projection of data with nadir " // &
1209	"location at:", lonNadir, latNadir
1210	ELSE
1211	msg = "Projection '" // TRIM(projN) // "' not ready" // CHAR(10) // " available ones: " // &
1212	"'ctsarea', 'lon/lat', 'nadir-sat'"
1213	CALL ErrMsg(msg,fname,-1)
1214	END IF
1215
1216	area = 0.
1217	sumRESx = 0.
1218	sumRESy = 0.
1219	meanwctr = 0.
1220	meanvnctr = 0.
1221	meanvxctr = 0.
1222	xcorr = 0.
1223	ycorr = 0.
1224
1225	nval = fillval64
1226	xval = -fillval64
1227
1228	ip = 1
1229	DO i=1,dx
1230	DO j=1,dy
1231	IF (poly(i,j)) THEN
1232	polypts(ip,1) = i
1233	polypts(ip,2) = j
1234	polyvals(ip) = values(i,j)
1235	SELECT CASE (TRIM(projS))
1236	CASE ('ctsarea')
1237	! Constant Area
1238	xcorr(i,j) = 1.
1239	ycorr(i,j) = 1.
1240	CASE ('lon/lat')
1241	! Area as fixed yres and sinus-lat varying for xres
1242	! IF (KIND(xcorr(i,j)) == KIND(1.d0)) THEN
1243	! xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
1244	! ELSE
1245	xcorr(i,j) = ABS(SIN(lon(i,j)*DegRad))
1246	! END IF
1247	ycorr(i,j) = 1.
1248	CASE ('nadir-sat')
1249	! Area from nadir resolution and degrading as we get far from satellite's nadir
1250	! GOES-E: 0 N, 75 W
1251	! IF (KIND(xcorr(i,j)) == KIND(1.d0)) THEN
1252	! xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
1253	! ELSE
1254	xcorr(i,j) = ABS(SIN(lon(i,j)*DegRad))
1255	! END IF
1256	ycorr(i,j) = 1.
1257	END SELECT
1258	area = area + xresxcorr(i,j)yres*ycorr(i,j)
1259	meanwctr(1) = meanwctr(1) + lon(i,j)xresxcorr(i,j)
1260	meanwctr(2) = meanwctr(2) + lat(i,j)yresycorr(i,j)
1261	IF (nval > values(i,j)) THEN
1262	nvcoord(1) = i
1263	nvcoord(2) = j
1264	nval = values(i,j)
1265	END IF
1266	IF (xval < values(i,j)) THEN
1267	xvcoord(1) = i
1268	xvcoord(2) = j
1269	xval = values(i,j)
1270	END IF
1271	ip = ip + 1
1272	END IF
1273	END DO
1274	END DO
1275
1276	IF (dbg) THEN
1277	PRINT *,' grid_coord lon lat value _______ '
1278	DO ip=1, Npolypts
1279	PRINT *, polypts(ip,:), ';', lon(polypts(ip,1),polypts(ip,2)), lat(polypts(ip,1),polypts(ip,2)),&
1280	':', polyvals(ip)
1281	END DO
1282	END IF
1283
1284	sumRESx = xres*SUM(xcorr)
1285	sumRESy = yres*SUM(ycorr)
1286
1287	xxtrm = (/ MINVAL(polypts(:,1)), MAXVAL(polypts(:,1)) /)
1288	yxtrm = (/ MINVAL(polypts(:,2)), MAXVAL(polypts(:,2)) /)
1289	meanctr = (/ SUM(polypts(:,1))/Npolypts, SUM(polypts(:,2))/Npolypts /)
1290	meanwctr = (/ meanwctr(1)/sumRESx, meanwctr(2)/sumRESy /)
1291
1292	IF (dbg) THEN
1293	PRINT *,' mean grid center: ', meanctr, ' weighted mean center: ', meanwctr
1294	END IF
1295
1296	! Statistics of the values within the polygon
1297	CALL StatsR_K(Npolypts, polyvals, nval, xval, mval, m2val, stdval)
1298
1299	IF (dbg) THEN
1300	PRINT *,' minimum value: ', nval, ' maximum:', xval, ' mean:', mval
1301	PRINT *,' coor. minimum: ', nvcoord
1302	PRINT *,' coor. maximum: ', xvcoord
1303	END IF
1304
1305	! Mean center weighted to minimum and maximum values
1306	! IF (KIND(polyvals(1)) == KIND(1.d0)) THEN
1307	! distvn = DABS(polyvals - nval)
1308	! distvx = DABS(polyvals - xval)
1309	! ELSE
1310	distvn = ABS(polyvals - nval)
1311	distvx = ABS(polyvals - xval)
1312	! END IF
1313
1314	meanvnctr = 0.
1315	meanvxctr = 0.
1316	sumDISTnlon = 0.
1317	sumDISTnlat = 0.
1318	sumDISTxlon = 0.
1319	sumDISTxlat = 0.
1320
1321	ip = 1
1322	DO i=1,dx
1323	DO j=1,dy
1324	IF (poly(i,j)) THEN
1325	meanvnctr(1) = meanvnctr(1)+lon(i,j)distvn(ip)xres*xcorr(i,j)
1326	meanvnctr(2) = meanvnctr(2)+lat(i,j)distvn(ip)yres*ycorr(i,j)
1327
1328	meanvxctr(1) = meanvxctr(1)+lon(i,j)distvx(ip)xres*xcorr(i,j)
1329	meanvxctr(2) = meanvxctr(2)+lat(i,j)distvx(ip)yres*ycorr(i,j)
1330
1331	sumDISTnlon = sumDISTnlon + distvn(ip)xresxcorr(i,j)
1332	sumDISTnlat = sumDISTnlat + distvn(ip)yresycorr(i,j)
1333	sumDISTxlon = sumDISTxlon + distvx(ip)xresxcorr(i,j)
1334	sumDISTxlat = sumDISTxlat + distvx(ip)yresycorr(i,j)
1335
1336	ip = ip + 1
1337	END IF
1338	END DO
1339	END DO
1340
1341	meanvnctr = (/ meanvnctr(1)/sumDISTnlon, meanvnctr(2)/sumDISTnlat /)
1342	meanvxctr = (/ meanvxctr(1)/sumDISTxlon, meanvxctr(2)/sumDISTxlat /)
1343
1344	IF (dbg) THEN
1345	PRINT *,' mean center to minimum: ', meanvnctr, ' to maximum: ', meanvxctr
1346	END IF
1347
1348	! Quantiles of the values within the polygon
1349	quant = -9999.d0
1350	IF (Npolypts > Nquant) THEN
1351	CALL quantilesR_K(Npolypts, polyvals, Nquant, quant)
1352	ELSE
1353	CALL SortR_K(polyvals, Npolypts)
1354	END IF
1355
1356	DEALLOCATE (polypts)
1357	DEALLOCATE (polyvals)
1358
1359	RETURN
1360
1361	END SUBROUTINE polygon_properties
1362
1363	SUBROUTINE polygons_t(dbg, dx, dy, dt, boolmatt, polys, Npoly)
1364	! Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!
1365
1366	IMPLICIT NONE
1367
1368	INTEGER, INTENT(in) :: dx, dy, dt
1369	LOGICAL, DIMENSION(dx,dy,dt), INTENT(in) :: boolmatt
1370	LOGICAL, INTENT(in) :: dbg
1371	INTEGER, DIMENSION(dt), INTENT(out) :: Npoly
1372	INTEGER, DIMENSION(dx,dy,dt), INTENT(out) :: polys
1373
1374	! Local
1375	INTEGER :: i,it
1376
1377	!!!!!!! Variables
1378	! dx,dy: spatial dimensions of the space
1379	! boolmatt: boolean matrix tolook for the polygons (.TRUE. based)
1380	! polys: found polygons
1381	! Npoly: number of polygons found
1382
1383	fname = 'polygons'
1384
1385	IF (dbg) PRINT *,TRIM(fname)
1386
1387	polys = -1
1388	Npoly = 0
1389
1390	DO it=1,dt
1391	IF (ANY(boolmatt(:,:,it))) THEN
1392	IF (dbg) THEN
1393	PRINT *,' it:', it, ' num. TRUE:', COUNT(boolmatt(:,:,it)), 'bool _______'
1394	DO i=1,dx
1395	PRINT *,boolmatt(i,:,it)
1396	END DO
1397	END IF
1398	CALL polygons(dbg, dx, dy, boolmatt(:,:,it), polys(:,:,it), Npoly(it))
1399	ELSE
1400	IF (dbg) THEN
1401	PRINT *,' it:', it, " without '.TRUE.' values skipiing it!!"
1402	END IF
1403	END IF
1404	END DO
1405
1406	END SUBROUTINE polygons_t
1407
1408	SUBROUTINE polygons(dbg, dx, dy, boolmat, polys, Npoly)
1409	! Subroutine to search the polygons of a boolean field. FORTRAN based. 1st = 1!
1410
1411	IMPLICIT NONE
1412
1413	INTEGER, INTENT(in) :: dx, dy
1414	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: boolmat
1415	LOGICAL, INTENT(in) :: dbg
1416	INTEGER, INTENT(out) :: Npoly
1417	INTEGER, DIMENSION(dx,dy), INTENT(out) :: polys
1418
1419	! Local
1420	INTEGER :: i, j, ip, ipp, Nppt
1421	INTEGER :: ierr
1422	INTEGER, DIMENSION(:,:), ALLOCATABLE :: borders
1423	LOGICAL, DIMENSION(dx,dy) :: isborder, isbordery
1424	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: paths
1425	INTEGER :: Npath
1426	INTEGER, DIMENSION(:), ALLOCATABLE :: Nptpaths
1427	INTEGER, DIMENSION(2) :: xtrx, xtry, meanpth
1428	INTEGER :: Nvertx, Npts
1429	INTEGER, DIMENSION(:,:), ALLOCATABLE :: vertxs, points
1430	LOGICAL, DIMENSION(:), ALLOCATABLE :: isin
1431
1432	!!!!!!! Variables
1433	! dx,dy: spatial dimensions of the space
1434	! boolmat: boolean matrix tolook for the polygons (.TRUE. based)
1435	! polys: found polygons
1436	! Npoly: number of polygons found
1437
1438	fname = 'polygons'
1439
1440	polys = -1
1441
1442	! The mathematical maximum woiuld be dx*dy/4, but let's be optimistic... (sorry Jero)
1443	Nppt = dx*dy/10
1444
1445	IF (ALLOCATED(borders)) DEALLOCATE(borders)
1446	ALLOCATE(borders(Nppt,2), STAT=ierr)
1447	msg = "Problems allocating matrix 'borders'"
1448	CALL ErrMsg(msg, fname, ierr)
1449
1450	IF (ALLOCATED(paths)) DEALLOCATE(paths)
1451	ALLOCATE(paths(Nppt,Nppt,2), STAT=ierr)
1452	msg = "Problems allocating matrix 'paths'"
1453	CALL ErrMsg(msg, fname, ierr)
1454
1455	IF (ALLOCATED(Nptpaths)) DEALLOCATE(Nptpaths)
1456	ALLOCATE(Nptpaths(Nppt), STAT=ierr)
1457	msg = "Problems allocating matrix 'Nptpaths'"
1458	CALL ErrMsg(msg, fname, ierr)
1459
1460	! Filling with the points of all the space with .TRUE.
1461	Npts = COUNT(boolmat)
1462
1463	IF (ALLOCATED(points)) DEALLOCATE(points)
1464	ALLOCATE(points(Npts,2), STAT=ierr)
1465	msg = "Problems allocating matrix 'points'"
1466	CALL ErrMsg(msg, fname, ierr)
1467
1468	! We only want to localize that points 'inside'
1469	ip = 1
1470	DO i=1, dx
1471	DO j=1, dy
1472	IF (boolmat(i,j)) THEN
1473	points(ip,1) = i
1474	points(ip,2) = j
1475	ip = ip + 1
1476	END IF
1477	END DO
1478	END DO
1479
1480	CALL borders_matrixL(dx, dy, Nppt, boolmat, borders, isborder, isbordery)
1481	CALL paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)
1482
1483	Npoly = Npath
1484
1485	DO ip=1, Npath
1486	IF (ALLOCATED(vertxs)) DEALLOCATE(vertxs)
1487	ALLOCATE(vertxs(Nptpaths(ip),2))
1488	msg = "Problems allocating matrix 'vertxs'"
1489	CALL ErrMsg(msg, fname, ierr)
1490
1491	IF (ALLOCATED(isin)) DEALLOCATE(isin)
1492	ALLOCATE(isin(Npts), STAT=ierr)
1493	msg = "Problems allocating matrix 'isin'"
1494	CALL ErrMsg(msg, fname, ierr)
1495
1496	isin = .FALSE.
1497
1498	IF (dbg) PRINT *, ' path:', ip, ' N pts:', Nptpaths(ip)
1499
1500	CALL path_properties(dx, dy, boolmat, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), xtrx, xtry, &
1501	meanpth, 'y', Nvertx, vertxs)
1502
1503	IF (dbg) THEN
1504	PRINT *, ' properties _______'
1505	PRINT *, ' x-extremes:', xtrx
1506	PRINT *, ' y-extremes:', xtry
1507	PRINT *, ' center mean:', meanpth
1508	PRINT *, ' y-vertexs:', Nvertx,' ________'
1509	DO i=1, Nvertx
1510	PRINT *,' ',i,':',vertxs(i,:)
1511	END DO
1512	END IF
1513
1514	CALL gridpoints_InsidePolygon(dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), Nvertx,&
1515	xtrx, xtry, vertxs, Npts, points, isin)
1516
1517	! Filling polygons
1518	DO ipp=1, Npts
1519	IF (isin(ipp)) polys(points(ipp,1),points(ipp,2)) = ip
1520	END DO
1521
1522	IF (dbg) THEN
1523	PRINT *,' boolmat isborder isbordery polygon (',xtrx(1),',',xtry(1),')x(',xtrx(2),',',xtry(2), &
1524	') _______'
1525	DO i=xtrx(1), xtrx(2)
1526	PRINT *,i,':',boolmat(i,xtry(1):xtry(2)), ' border ', isborder(i,xtry(1):xtry(2)), &
1527	' isbordery ', isbordery(i,xtry(1):xtry(2)), ' polygon ', polys(i,xtry(1):xtry(2))
1528	END DO
1529	END IF
1530
1531	END DO
1532
1533	! Cleaning polygons matrix of not-used and paths around holes
1534	CALL clean_polygons(dx, dy, boolmat, polys, Npoly, dbg)
1535
1536	DEALLOCATE (borders)
1537	DEALLOCATE (Nptpaths)
1538	DEALLOCATE (paths)
1539	DEALLOCATE (vertxs)
1540	DEALLOCATE (points)
1541	DEALLOCATE (isin)
1542
1543	RETURN
1544
1545	END SUBROUTINE polygons
1546
1547	SUBROUTINE clean_polygons(dx, dy, Lmat, pols, Npols, dbg)
1548	! Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
1549
1550	IMPLICIT NONE
1551
1552	INTEGER, INTENT(in) :: dx, dy
1553	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat
1554	INTEGER, INTENT(inout) :: Npols
1555	INTEGER, DIMENSION(dx,dy), INTENT(inout) :: pols
1556	LOGICAL, INTENT(in) :: dbg
1557
1558	! Local
1559	INTEGER :: i,j,ip,iprm
1560	INTEGER, DIMENSION(Npols) :: origPol, NotPol, neigPol
1561	INTEGER :: ispol, NnotPol
1562	CHARACTER(len=4) :: ISa
1563
1564	!!!!!!! Variables
1565	! dx, dy: size of the space
1566	! Lmat: original bolean matrix from which the polygons come from
1567	! Npols: original number of polygons
1568	! pols: polygons space
1569
1570	fname = 'clean_polygons'
1571	IF (dbg) PRINT *," At '" // TRIM(fname) // "' ..."
1572
1573	origPol = -1
1574
1575	! Looking for polygons already in space
1576	NnotPol = 0
1577	DO ip=1, Npols
1578	ispol = COUNT(pols-ip == 0)
1579	IF (ispol > 0) THEN
1580	origPol(ip) = ip
1581	ELSE
1582	NnotPol = NnotPol + 1
1583	NotPol(NnotPol) = ip
1584	neigPol(NnotPol) = -1
1585	END IF
1586	END DO
1587
1588	IF (dbg) THEN
1589	PRINT *,' It should be:', Npols, ' polygons, but already there are:', Npols - NnotPol
1590	PRINT *,' Polygons to remove:', NotPol(1:NnotPol)
1591	END IF
1592
1593	! Looking for the hole border of a polygon. This is identify as such polygon point which along
1594	! y-axis has NpolygonA, Npolygon, .FALSE.
1595	DO i=1,dx
1596	DO j=2,dy-1
1597	IF ( (pols(i,j-1) /= pols(i,j) .AND. pols(i,j+1) == -1) .AND. (COUNT(NotPol-pols(i,j)==0)==0) &
1598	.AND. (pols(i,j) /= -1) .AND. (pols(i,j-1) /= -1)) THEN
1599	IF (dbg) PRINT *,' Polygon:', pols(i,j), ' to be removed at point (',i,',',j,'); j-1:', &
1600	pols(i,j-1), ' j:', pols(i,j), ' j+1:', pols(i,j+1)
1601	NnotPol = NnotPol + 1
1602	NotPol(NnotPol) = pols(i,j)
1603	neigPol(NnotPol) = pols(i,j-1)
1604	END IF
1605	END DO
1606	END DO
1607
1608	IF (dbg) THEN
1609	PRINT *,' It should be:', Npols, ' polygons, but already there are:', Npols - NnotPol
1610	PRINT *,' Polygons to remove after looking for fake border-of-hole polygons _______'
1611	DO i=1, NnotPol
1612	PRINT *, ' Polygon:', NotPol(i), ' to be replaced by:', neigPol(i)
1613	END DO
1614	END IF
1615
1616	! Removing polygons
1617	DO iprm=1, NnotPol
1618	IF (neigPol(iprm) == -1) THEN
1619	WHERE (pols == NotPol(iprm))
1620	pols = -1
1621	END WHERE
1622	IF (dbg) THEN
1623	PRINT *,' removing polygon:', NotPol(iprm)
1624	END IF
1625	ELSE
1626	WHERE (pols == NotPol(iprm))
1627	pols = neigPol(iprm)
1628	END WHERE
1629	IF (dbg) THEN
1630	PRINT *,' replacing polygon:', NotPol(iprm), ' by:', neigPol(iprm)
1631	END IF
1632	END IF
1633	END DO
1634
1635	! Re-numbering (descending values)
1636	DO i = 1, NnotPol
1637	iprm = MAXVAL(NotPol(1:NnotPol))
1638	WHERE(pols > iprm)
1639	pols = pols - 1
1640	END WHERE
1641	j = Index1DArrayI(NotPol, NnotPol, iprm)
1642	NotPol(j) = -9
1643	END DO
1644
1645	Npols = Npols - NnotPol
1646
1647	RETURN
1648
1649	END SUBROUTINE clean_polygons
1650
1651	SUBROUTINE path_properties(dx, dy, Lmat, Nptspth, pth, xxtrm, yxtrm, meanctr, axs, Nvrtx, vrtxs)
1652	! Subroutine to determine the properties of a path:
1653	! extremes: minimum and maximum of the path along x,y axes
1654	! meancenter: center from the mean of the coordinates of the paths locations
1655	! vertexs: path point, without neighbours along a given axis
1656
1657	IMPLICIT NONE
1658
1659	INTEGER, INTENT(in) :: dx, dy, Nptspth
1660	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat
1661	INTEGER, DIMENSION(Nptspth,2), INTENT(in) :: pth
1662	CHARACTER, INTENT(in) :: axs
1663	INTEGER, DIMENSION(2), INTENT(out) :: meanctr, xxtrm, yxtrm
1664	INTEGER, INTENT(out) :: Nvrtx
1665	INTEGER, DIMENSION(Nptspth,2), INTENT(out) :: vrtxs
1666
1667	! Local
1668	INTEGER :: i, ip, jp
1669	INTEGER :: neig1, neig2
1670
1671	!!!!!!! Variables
1672	! dx,dy: size of the space
1673	! Lmat: original matrix of logical values for the path
1674	! Nptspth: number of points of the path
1675	! pth: path coordinates (clockwise)
1676	! axs: axis of finding the vertex
1677	! [x/y]xtrm: minimum and maximum coordinates of the path
1678	! meanctr: center from the mean of the coordinates of the path
1679	! Nvrtx: Number of vertexs of the path
1680	! vrtxs: coordinates of the vertexs
1681
1682	fname = 'path_properties'
1683
1684	vrtxs = -1
1685	Nvrtx = 0
1686
1687	xxtrm = (/ MINVAL(pth(:,1)), MAXVAL(pth(:,1)) /)
1688	yxtrm = (/ MINVAL(pth(:,2)), MAXVAL(pth(:,2)) /)
1689	meanctr = (/ SUM(pth(:,1))/Nptspth, SUM(pth(:,2))/Nptspth /)
1690
1691	IF (axs == 'x' .OR. axs == 'X') THEN
1692	! Looking vertexs along x-axis
1693	DO i=1, Nptspth
1694	ip = pth(i,1)
1695	jp = pth(i,2)
1696	neig1 = 0
1697	neig2 = 0
1698	! W-point
1699	IF (ip == 1) THEN
1700	neig1 = -1
1701	ELSE
1702	IF (.NOT.Lmat(ip-1,jp)) neig1 = -1
1703	END IF
1704	! E-point
1705	IF (ip == dx) THEN
1706	neig2 = -1
1707	ELSE
1708	IF (.NOT.Lmat(ip+1,jp)) neig2 = -1
1709	END IF
1710
1711	IF (neig1 == -1 .AND. neig2 == -1) THEN
1712	Nvrtx = Nvrtx + 1
1713	vrtxs(Nvrtx,:) = (/ip,jp/)
1714	END IF
1715	END DO
1716	ELSE IF (axs == 'y' .OR. axs == 'Y') THEN
1717	! Looking vertexs along x-axis
1718	DO i=1, Nptspth
1719	ip = pth(i,1)
1720	jp = pth(i,2)
1721
1722	neig1 = 0
1723	neig2 = 0
1724	! S-point
1725	IF (jp == 1) THEN
1726	neig1 = -1
1727	ELSE
1728	IF (.NOT.Lmat(ip,jp-1)) neig1 = -1
1729	END IF
1730	! N-point
1731	IF (jp == dy) THEN
1732	neig2 = -1
1733	ELSE
1734	IF (.NOT.Lmat(ip,jp+1)) neig2 = -1
1735	END IF
1736
1737	IF (neig1 == -1 .AND. neig2 == -1) THEN
1738	Nvrtx = Nvrtx + 1
1739	vrtxs(Nvrtx,:) = (/ ip, jp /)
1740	END IF
1741	END DO
1742	ELSE
1743	msg = "Axis '" // axs // "' not available" // CHAR(10) // " Available ones: 'x', 'X', 'y, 'Y'"
1744	CALL ErrMsg(msg, fname, -1)
1745	END IF
1746
1747	RETURN
1748
1749	END SUBROUTINE path_properties
1750
1751	SUBROUTINE gridpoints_InsidePolygon(dx, dy, isbrdr, Npath, path, Nvrtx, xpathxtrm, ypathxtrm, &
1752	vrtxs, Npts, pts, inside)
1753	! Subroutine to determine if a series of grid points are inside a polygon following ray casting algorithm
1754	! FROM: https://en.wikipedia.org/wiki/Point_in_polygon
1755
1756	IMPLICIT NONE
1757
1758	INTEGER, INTENT(in) :: dx,dy,Npath,Nvrtx,Npts
1759	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isbrdr
1760	INTEGER, DIMENSION(Npath,2), INTENT(in) :: path
1761	INTEGER, DIMENSION(2), INTENT(in) :: xpathxtrm, ypathxtrm
1762	INTEGER, DIMENSION(Npath,2) :: vrtxs
1763	INTEGER, DIMENSION(Npts,2), INTENT(in) :: pts
1764	LOGICAL, DIMENSION(Npts), INTENT(out) :: inside
1765
1766	! Local
1767	INTEGER :: i,j,ip,ix,iy
1768	INTEGER :: Nintersecs, isvertex, ispath
1769	INTEGER :: ierr
1770	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: halo_brdr
1771	INTEGER :: Nbrbrdr
1772
1773	!!!!!!! Variables
1774	! dx,dy: space size
1775	! Npath: number of points of the path of the polygon
1776	! path: path of the polygon
1777	! isbrdr: boolean matrix of the space wqith .T. on polygon border
1778	! Nvrtx: number of vertexs of the path
1779	! [x/y]pathxtrm extremes of the path
1780	! vrtxs: vertexs of the path along y-axis
1781	! Npts: number of points
1782	! pts: points to look for
1783	! inside: vector wether point is inside or not (coincident to a border is inside)
1784
1785	fname = 'gridpoints_InsidePolygon'
1786
1787	! Creation of a 1-grid point larger matrix to deal with points reaching the limits
1788	IF (ALLOCATED(halo_brdr)) DEALLOCATE(halo_brdr)
1789	ALLOCATE(halo_brdr(dx+2,dy+2), STAT=ierr)
1790	msg = "Problems allocating matrix 'halo_brdr'"
1791	CALL ErrMsg(msg, fname, ierr)
1792	halo_brdr = .FALSE.
1793
1794	DO i=1,dx
1795	halo_brdr(i+1,2:dy+1) = isbrdr(i,:)
1796	END DO
1797
1798	inside = .FALSE.
1799
1800	DO ip=1,Npts
1801	Nintersecs = 0
1802	ix = pts(ip,1)
1803	iy = pts(ip,2)
1804	! Point might be outside path range...
1805	IF (ix >= xpathxtrm(1) .AND. ix <= xpathxtrm(2) .AND. iy >= ypathxtrm(1) .AND. &
1806	iy <= ypathxtrm(2)) THEN
1807
1808	! It is a border point?
1809	ispath = index_list_coordsI(Npath, path, (/ix,iy/))
1810	IF (isbrdr(ix,iy) .AND. (ispath /= -1)) THEN
1811	inside(ip) = .TRUE.
1812	CYCLE
1813	END IF
1814
1815	! Looking along y-axis
1816	! Accounting for consecutives borders
1817	Nbrbrdr = 0
1818	DO j=MAX(1,ypathxtrm(1)-1),iy-1
1819	! Only counting that borders that are not vertexs
1820	ispath = index_list_coordsI(Npath, path, (/ix,j/))
1821	isvertex = index_list_coordsI(Npath, vrtxs, (/ix,j/))
1822
1823	IF (halo_brdr(ix+1,j+1) .AND. (ispath /= -1) .AND. (isvertex == -1) ) Nintersecs = Nintersecs + 1
1824	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
1825	Nbrbrdr = Nbrbrdr + 1
1826	ELSE
1827	! Will remove that consecutive borders above 2
1828	IF (Nbrbrdr /= 0) THEN
1829	Nintersecs = Nintersecs - MAX(Nbrbrdr-1, 0)
1830	Nbrbrdr = 0
1831	END IF
1832	END IF
1833	END DO
1834	IF (MOD(Nintersecs,2) /= 0) inside(ip) = .TRUE.
1835	END IF
1836
1837	END DO
1838
1839	RETURN
1840
1841	END SUBROUTINE gridpoints_InsidePolygon
1842
1843	SUBROUTINE look_clockwise_borders(dx,dy,Nbrdrs,brdrs,gbrdr,isbrdr,ix,iy,dbg,xf,yf,iff)
1844	! Subroutine to look clock-wise for a next point within a collection of borders (limits of a region)
1845
1846	IMPLICIT NONE
1847
1848	INTEGER, INTENT(in) :: dx, dy, Nbrdrs, ix, iy
1849	INTEGER, DIMENSION(Nbrdrs,2), INTENT(in) :: brdrs
1850	LOGICAL, DIMENSION(Nbrdrs), INTENT(in) :: gbrdr
1851	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isbrdr
1852	LOGICAL, INTENT(in) :: dbg
1853	INTEGER, INTENT(out) :: xf, yf, iff
1854
1855	! Local
1856	INTEGER :: isch
1857	CHARACTER(len=2), DIMENSION(8) :: Lclock
1858	INTEGER, DIMENSION(8,2) :: spt
1859	INTEGER :: iif, jjf
1860
1861	!!!!!!! Variables
1862	! dx, dy: 2D shape ot the space
1863	! Nbrdrs: number of brdrs found in this 2D space
1864	! brdrs: list of coordinates of the borders
1865	! gbrdr: accounts for the use if the given border point
1866	! isbrdr: accounts for the matrix of the point is a border or not
1867	! ix,iy: coordinates of the point to start to find for
1868	! xf,yf: coordinates of the found point
1869	! iff: position of the border found within the list of borders
1870
1871	fname = 'look_clockwise_borders'
1872
1873	! Looking clock-wise assuming that one starts from the westernmost point
1874
1875	! Label of the search
1876	lclock = (/ 'W ', 'NW', 'N ', 'NE', 'E ', 'SE', 'S ', 'SW' /)
1877	! Transformation to apply
1878	!spt = (/ (/-1,0/), (/-1,1/), (/0,1/), (/1,1/), (/1,0/), (/1,-1/), (/0,-1/), (/-1,-1/) /)
1879	spt(:,1) = (/ -1, -1, 0, 1, 1, 1, 0, -1 /)
1880	spt(:,2) = (/ 0, 1, 1, 1, 0, -1, -1, -1 /)
1881
1882	xf = -1
1883	yf = -1
1884	DO isch=1, 8
1885	! clock-wise search
1886	IF (spt(isch,1) >= 0) THEN
1887	iif = MIN(dx,ix+spt(isch,1))
1888	ELSE
1889	iif = MAX(1,ix+spt(isch,1))
1890	END IF
1891	IF (spt(isch,2) >= 0) THEN
1892	jjf = MIN(dy,iy+spt(isch,2))
1893	ELSE
1894	jjf = MAX(1,iy+spt(isch,2))
1895	END IF
1896	iff = index_list_coordsI(Nbrdrs, brdrs,(/iif,jjf/))
1897	IF (iff > 0) THEN
1898	IF (dbg) PRINT *,' ' // lclock(isch) // '-point:', iif,jjf, ':', iff, 'is',isbrdr(iif,jjf), &
1899	'got',gbrdr(iff)
1900	IF (isbrdr(iif,jjf) .AND. .NOT.gbrdr(iff)) THEN
1901	xf = iif
1902	yf = jjf
1903	EXIT
1904	END IF
1905	END IF
1906	END DO
1907
1908	RETURN
1909
1910	END SUBROUTINE look_clockwise_borders
1911
1912	SUBROUTINE borders_matrixL(dx,dy,dxy,Lmat,brdrs,isbrdr,isbrdry)
1913	! Subroutine to provide the borders of a logical array (interested in .TRUE.)
1914
1915	IMPLICIT NONE
1916
1917	INTEGER, INTENT(in) :: dx,dy,dxy
1918	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat
1919	INTEGER, DIMENSION(dxy,2), INTENT(out) :: brdrs
1920	LOGICAL, DIMENSION(dx,dy), INTENT(out) :: isbrdr, isbrdry
1921
1922	! Local
1923	INTEGER :: i,j,ib
1924
1925	!!!!!!! Variables
1926	! dx,dy: size of the space
1927	! dxy: maximum number of border points
1928	! Lmat: Matrix to look for the borders
1929	! brdrs: list of coordinates of the borders
1930	! isbrdr: matrix with .T./.F. wether the given matrix point is a border or not
1931	! isbrdry: matrix with .T./.F. wether the given matrix point is a border or not only along y-axis
1932
1933	fname = 'borders_matrixL'
1934
1935	isbrdr = .FALSE.
1936	brdrs = -1
1937	ib = 1
1938
1939	! Starting with the borders. If a given point is TRUE it is a path-vertex
1940	! Along y-axis
1941	DO i=1, dx
1942	IF (Lmat(i,1) .AND. .NOT.isbrdr(i,1)) THEN
1943	brdrs(ib,1) = i
1944	brdrs(ib,2) = 1
1945	isbrdr(i,1) = .TRUE.
1946	ib=ib+1
1947	END IF
1948	IF (Lmat(i,dy) .AND. .NOT.isbrdr(i,dy)) THEN
1949	brdrs(ib,1) = i
1950	brdrs(ib,2) = dy
1951	isbrdr(i,dy) = .TRUE.
1952	ib=ib+1
1953	END IF
1954	END DO
1955	! Along x-axis
1956	DO j=1, dy
1957	IF (Lmat(1,j) .AND. .NOT.isbrdr(1,j)) THEN
1958	brdrs(ib,1) = 1
1959	brdrs(ib,2) = j
1960	isbrdr(1,j) = .TRUE.
1961	ib=ib+1
1962	END IF
1963	IF (Lmat(dx,j) .AND. .NOT.isbrdr(dx,j)) THEN
1964	brdrs(ib,1) = dx
1965	brdrs(ib,2) = j
1966	isbrdr(dx,j) = .TRUE.
1967	ib=ib+1
1968	END IF
1969	END DO
1970
1971	isbrdry = isbrdr
1972
1973	! Border as that when looking on x-axis points with Lmat(i) /= Lmat(i+1)
1974	DO i=1, dx-1
1975	DO j=1, dy-1
1976	IF ( Lmat(i,j) .NEQV. Lmat(i+1,j) ) THEN
1977	IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN
1978	brdrs(ib,1) = i
1979	brdrs(ib,2) = j
1980	isbrdr(i,j) = .TRUE.
1981	ib=ib+1
1982	ELSE IF (Lmat(i+1,j) .AND. .NOT.isbrdr(i+1,j)) THEN
1983	brdrs(ib,1) = i+1
1984	brdrs(ib,2) = j
1985	isbrdr(i+1,j) = .TRUE.
1986	ib=ib+1
1987	END IF
1988	END IF
1989	! y-axis
1990	IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN
1991	IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN
1992	brdrs(ib,1) = i
1993	brdrs(ib,2) = j
1994	isbrdr(i,j) = .TRUE.
1995	isbrdry(i,j) = .TRUE.
1996	ib=ib+1
1997	ELSE IF (Lmat(i,j+1) .AND. .NOT.isbrdr(i,j+1)) THEN
1998	brdrs(ib,1) = i
1999	brdrs(ib,2) = j+1
2000	isbrdr(i,j+1) = .TRUE.
2001	isbrdry(i,j+1) = .TRUE.
2002	ib=ib+1
2003	END IF
2004	END IF
2005	END DO
2006	END DO
2007
2008	DO i=1, dx-1
2009	DO j=1, dy-1
2010	! y-axis
2011	IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN
2012	IF (Lmat(i,j)) THEN
2013	isbrdry(i,j) = .TRUE.
2014	ELSE IF (Lmat(i,j+1)) THEN
2015	isbrdry(i,j+1) = .TRUE.
2016	END IF
2017	END IF
2018	END DO
2019	END DO
2020	! only y-axis adding bands of 2 grid points
2021	DO i=1, dx-1
2022	DO j=2, dy-2
2023	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
2024	IF (Lmat(i,j)) THEN
2025	isbrdry(i,j) = .TRUE.
2026	isbrdry(i,j+1) = .TRUE.
2027	END IF
2028	END IF
2029	END DO
2030	END DO
2031
2032	RETURN
2033
2034	END SUBROUTINE borders_matrixL
2035
2036	SUBROUTINE paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)
2037	! Subroutine to search the paths of a border field.
2038
2039	IMPLICIT NONE
2040
2041	INTEGER, INTENT(in) :: dx, dy, Nppt
2042	LOGICAL, INTENT(in) :: dbg
2043	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isborder
2044	INTEGER, DIMENSION(Nppt,2), INTENT(in) :: borders
2045	INTEGER, DIMENSION(Nppt,Nppt,2), INTENT(out) :: paths
2046	INTEGER, INTENT(out) :: Npath
2047	INTEGER, DIMENSION(Nppt), INTENT(out) :: Nptpaths
2048
2049	! Local
2050	INTEGER :: i,j,k,ib
2051	INTEGER :: ierr
2052	INTEGER :: Nbrdr
2053	LOGICAL, DIMENSION(:), ALLOCATABLE :: gotbrdr, emptygotbrdr
2054	INTEGER :: iipth, ipath, ip, Nptspath
2055	INTEGER :: iib, jjb, iip, ijp, iif, jjf, iff
2056	LOGICAL :: found, finishedstep
2057
2058	!!!!!!! Variables
2059	! dx,dy: spatial dimensions of the space
2060	! Nppt: possible number of paths and points that the paths can have
2061	! isborder: boolean matrix which provide the borders of the polygon
2062	! borders: coordinates of the borders of the polygon
2063	! paths: coordinates of each found path
2064	! Npath: number of paths found
2065	! Nptpaths: number of points per path
2066
2067	fname = 'paths_border'
2068
2069	! Sarting matrix
2070	paths = -1
2071	Npath = 0
2072	Nptspath = 0
2073	Nptpaths = -1
2074
2075	ib=1
2076	finishedstep = .FALSE.
2077
2078	! Number of border points
2079	DO ib=1, Nppt
2080	IF (borders(ib,1) == -1 ) EXIT
2081	END DO
2082	Nbrdr = ib-1
2083
2084	IF (dbg) THEN
2085	PRINT *,' borders _______'
2086	DO i=1,Nbrdr
2087	PRINT *,' ',i,':',borders(i,:)
2088	END DO
2089	END IF
2090
2091	! Matrix which keeps track if a border point has been located
2092	IF (ALLOCATED(gotbrdr)) DEALLOCATE(gotbrdr)
2093	ALLOCATE(gotbrdr(Nbrdr), STAT=ierr)
2094	msg = "Problems allocating matrix 'gotbrdr'"
2095	CALL ErrMsg(msg, fname, ierr)
2096	IF (ALLOCATED(emptygotbrdr)) DEALLOCATE(emptygotbrdr)
2097	ALLOCATE(emptygotbrdr(Nbrdr), STAT=ierr)
2098	msg = "Problems allocating matrix 'emptygotbrdr'"
2099	CALL ErrMsg(msg, fname, ierr)
2100
2101	gotbrdr = .FALSE.
2102	emptygotbrdr = .FALSE.
2103
2104	! Starting the fun...
2105
2106	! Looking along the lines and when a border is found, starting from there in a clock-wise way
2107	iipth = 1
2108	ipath = 1
2109	DO ib=1,Nbrdr
2110	iib = borders(iipth,1)
2111	jjb = borders(iipth,2)
2112	! Starting new path
2113	newpath: IF (.NOT.gotbrdr(iipth)) THEN
2114	ip = 1
2115	Nptspath = 1
2116	paths(ipath,ip,:) = borders(iipth,:)
2117	gotbrdr(iipth) = .TRUE.
2118	! Looking for following clock-wise search
2119	! Not looking for W, because search starts from the W
2120	iip = iib
2121	ijp = jjb
2122	DO k=1,Nbrdr
2123	IF (dbg) PRINT *,ipath,'iip jip:', iip, ijp
2124	found = .FALSE.
2125	CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE.,iif,jjf,iff)
2126	IF (iif /= -1) THEN
2127	ip=ip+1
2128	paths(ipath,ip,:) = (/ iif,jjf /)
2129	found = .TRUE.
2130	gotbrdr(iff) = .TRUE.
2131	iip = iif
2132	ijp = jjf
2133	Nptspath = Nptspath + 1
2134	END IF
2135
2136	IF (dbg) THEN
2137	PRINT *,iib,jjb,' end of this round path:', ipath, '_____', gotbrdr
2138	DO i=1, Nptspath
2139	PRINT *,' ',i,':',paths(ipath,i,:)
2140	END DO
2141	END IF
2142	! If it is not found a next point, might be because it is a non-polygon related value
2143	IF (.NOT.found) THEN
2144	IF (dbg) PRINT *,'NOT FOUND !!!', gotbrdr
2145	! Are still there available borders?
2146	IF (ALL(gotbrdr) .EQV. .TRUE.) THEN
2147	finishedstep = .TRUE.
2148	Npath = ipath
2149	Nptpaths(ipath) = Nptspath
2150	EXIT
2151	ELSE
2152	Nptpaths(ipath) = Nptspath
2153	! Let's have a look if the previous points in the path have already some 'non-located' neighbourgs
2154	DO i=Nptspath,1,-1
2155	iip = paths(ipath,i,1)
2156	ijp = paths(ipath,i,2)
2157	CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE., iif, &
2158	jjf,iff)
2159	IF (iif /= -1 .AND. iff /= -1) THEN
2160	IF (dbg) PRINT *,' re-take path from point:', iif,',',jjf,' n-path:', iff
2161	found = .TRUE.
2162	iipth = index_list_coordsI(Nppt, borders, (/iip,ijp/))
2163	EXIT
2164	END IF
2165	END DO
2166	IF (.NOT.found) THEN
2167	! Looking for the next available border point for the new path
2168	DO i=1,Nbrdr
2169	IF (.NOT.gotbrdr(i)) THEN
2170	iipth = i
2171	EXIT
2172	END IF
2173	END DO
2174	IF (dbg) PRINT *,' Looking for next path starting at:', iipth, ' point:', &
2175	borders(iipth,:)
2176	ipath=ipath+1
2177	EXIT
2178	END IF
2179	END IF
2180	ELSE
2181	IF (dbg) PRINT *,' looking for next point...'
2182	END IF
2183	IF (finishedstep) EXIT
2184	END DO
2185	END IF newpath
2186	END DO
2187	Npath = ipath
2188	Nptpaths(ipath) = Nptspath
2189
2190	DEALLOCATE (gotbrdr)
2191	DEALLOCATE (emptygotbrdr)
2192
2193	RETURN
2194
2195	END SUBROUTINE paths_border
2196
2197	SUBROUTINE rand_sample(Nvals, Nsample, sample)
2198	! Subroutine to randomly sample a range of indices
2199
2200	IMPLICIT NONE
2201
2202	INTEGER, INTENT(in) :: Nvals, Nsample
2203	INTEGER, DIMENSION(Nsample), INTENT(out) :: sample
2204
2205	! Local
2206	INTEGER :: i, ind, jmax
2207	REAL, DIMENSION(Nsample) :: randv
2208	CHARACTER(len=50) :: fname
2209	LOGICAL :: found
2210	LOGICAL, DIMENSION(Nvals) :: issampled
2211	CHARACTER(len=256) :: msg
2212	CHARACTER(len=10) :: IS1, IS2
2213
2214	!!!!!!! Variables
2215	! Nvals: number of values
2216	! Nsamples: number of samples
2217	! sample: samnple
2218	fname = 'rand_sample'
2219
2220	IF (Nsample > Nvals) THEN
2221	WRITE(IS1,'(I10)')Nvals
2222	WRITE(IS2,'(I10)')Nsample
2223	msg = 'Sampling of ' // TRIM(IS1) // ' is too big for ' // TRIM(IS1) // 'values'
2224	CALL ErrMsg(msg, fname, -1)
2225	END IF
2226
2227	! Generation of random numbers always the same series during the whole program!
2228	CALL RANDOM_NUMBER(randv)
2229
2230	! Making sure that we do not repeat any value
2231	issampled = .FALSE.
2232
2233	DO i=1, Nsample
2234	! Generation of the index from the random numbers
2235	ind = MAX(INT(randv(i)*Nvals), 1)
2236
2237	IF (.NOT.issampled(ind)) THEN
2238	sample(i) = ind
2239	issampled(ind) = .TRUE.
2240	ELSE
2241	! Looking around the given index
2242	!PRINT *,' Index :', ind, ' already sampled!', issampled(ind)
2243	found = .FALSE.
2244	DO jmax=1, Nvals
2245	ind = MIN(ind+jmax, Nvals)
2246	IF (.NOT.issampled(ind)) THEN
2247	sample(i) = ind
2248	issampled(ind) = .TRUE.
2249	found = .TRUE.
2250	EXIT
2251	END IF
2252	ind = MAX(1, ind-jmax)
2253	IF (.NOT.issampled(ind)) THEN
2254	sample(i) = ind
2255	issampled(ind) = .TRUE.
2256	found = .TRUE.
2257	EXIT
2258	END IF
2259	END DO
2260	IF (.NOT.found) THEN
2261	msg = 'sampling could not be finished due to absence of available value!!'
2262	CALL ErrMsg(msg, fname, -1)
2263	END IF
2264	END IF
2265
2266	END DO
2267
2268	RETURN
2269
2270	END SUBROUTINE rand_sample
2271
2272	SUBROUTINE PrintQuantilesR_K(Nvals, vals, Nquants, qtvs, bspc)
2273	! Subroutine to print the quantiles of values REAL(r_k)
2274
2275	IMPLICIT NONE
2276
2277	INTEGER, INTENT(in) :: Nvals, Nquants
2278	REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals
2279	REAL(r_k), DIMENSION(Nquants), INTENT(in) :: qtvs
2280	CHARACTER(len=1000), OPTIONAL :: bspc
2281
2282	! Local
2283	INTEGER :: iq
2284	LOGICAL, DIMENSION(Nvals) :: search1, search2, search
2285	CHARACTER(len=6) :: RS1
2286	CHARACTER(len=50) :: fname
2287	CHARACTER(len=1000) :: bspcS
2288
2289	!!!!!!! Variables
2290	! vals: series of values
2291	! qtvs: values of the quantiles
2292	! bspc: base quantity of spaces
2293
2294	fname = 'PrintQuantilesR_K'
2295
2296	IF (PRESENT(bspc)) THEN
2297	bspcS = bspc
2298	ELSE
2299	bspcS = ' '
2300	END IF
2301
2302	DO iq=1, Nquants-1
2303
2304	WHERE (vals >= qtvs(iq))
2305	search1 = .TRUE.
2306	ELSEWHERE
2307	search1 = .FALSE.
2308	END WHERE
2309
2310	WHERE (vals < qtvs(iq+1))
2311	search2 = .TRUE.
2312	ELSEWHERE
2313	search2 = .FALSE.
2314	END WHERE
2315
2316	WHERE (search1 .AND. search2)
2317	search = .TRUE.
2318	ELSEWHERE
2319	search = .FALSE.
2320	END WHERE
2321
2322	WRITE(RS1, '(F6.2)')(iq)*100./(Nquants-1)
2323	PRINT *, TRIM(bspcS) // '[',iq,']', TRIM(RS1) // ' %:', qtvs(iq), 'N:', COUNT(search)
2324
2325	END DO
2326
2327	RETURN
2328
2329	END SUBROUTINE PrintQuantilesR_K
2330
2331	INTEGER FUNCTION FindMinimumR_K(x, dsize, Startv, Endv)
2332	! Function returns the location of the minimum in the section between Start and End.
2333
2334	IMPLICIT NONE
2335
2336	INTEGER, INTENT(in) :: dsize
2337	REAL(r_k), DIMENSION(dsize), INTENT(in) :: x
2338	INTEGER, INTENT(in) :: Startv, Endv
2339
2340	! Local
2341	REAL(r_k) :: Minimum
2342	INTEGER :: Location
2343	INTEGER :: i
2344
2345	Minimum = x(Startv) ! assume the first is the min
2346	Location = Startv ! record its position
2347	DO i = Startv+1, Endv ! start with next elements
2348	IF (x(i) < Minimum) THEN ! if x(i) less than the min?
2349	Minimum = x(i) ! Yes, a new minimum found
2350	Location = i ! record its position
2351	END IF
2352	END DO
2353
2354	FindMinimumR_K = Location ! return the position
2355
2356	END FUNCTION FindMinimumR_K
2357
2358	SUBROUTINE SwapR_K(a, b)
2359	! Subroutine swaps the values of its two formal arguments.
2360
2361	IMPLICIT NONE
2362
2363	REAL(r_k), INTENT(INOUT) :: a, b
2364	! Local
2365	REAL(r_k) :: Temp
2366
2367	Temp = a
2368	a = b
2369	b = Temp
2370
2371	END SUBROUTINE SwapR_K
2372
2373	SUBROUTINE SortR_K(x, Nx)
2374	! Subroutine receives an array x() r_K and sorts it into ascending order.
2375
2376	IMPLICIT NONE
2377
2378	INTEGER, INTENT(IN) :: Nx
2379	REAL(r_k), DIMENSION(Nx), INTENT(INOUT) :: x
2380
2381	! Local
2382	INTEGER :: i
2383	INTEGER :: Location
2384
2385	DO i = 1, Nx-1 ! except for the last
2386	Location = FindMinimumR_K(x, Nx-i+1, i, Nx) ! find min from this to last
2387	CALL SwapR_K(x(i), x(Location)) ! swap this and the minimum
2388	END DO
2389
2390	END SUBROUTINE SortR_K
2391
2392	SUBROUTINE quantilesR_K(Nvals, vals, Nquants, quants)
2393	! Subroutine to provide the quantiles of a given set of values of type real 'r_k'
2394
2395	IMPLICIT NONE
2396
2397	INTEGER, INTENT(in) :: Nvals, Nquants
2398	REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals
2399	REAL(r_k), DIMENSION(Nquants), INTENT(out) :: quants
2400
2401	! Local
2402	INTEGER :: i
2403	REAL(r_k) :: minv, maxv
2404	REAL(r_k), DIMENSION(Nvals) :: sortedvals
2405
2406	!!!!!!! Variables
2407	! Nvals: number of values
2408	! Rk: kind of real of the values
2409	! vals: values
2410	! Nquants: number of quants
2411	! quants: values at which the quantile start
2412
2413	minv = MINVAL(vals)
2414	maxv = MAXVAL(vals)
2415
2416	sortedvals = vals
2417	! Using from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
2418	CALL SortR_K(sortedvals, Nvals)
2419
2420	quants(1) = minv
2421	DO i=2, Nquants
2422	quants(i) = sortedvals(INT((i-1)*Nvals/Nquants))
2423	END DO
2424
2425	END SUBROUTINE quantilesR_K
2426
2427
2428	SUBROUTINE StatsR_K(Nvals, vals, minv, maxv, mean, mean2, stdev)
2429	! Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a
2430	! series of r_k numbers
2431
2432	IMPLICIT NONE
2433
2434	INTEGER, INTENT(in) :: Nvals
2435	REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals
2436	REAL(r_k), INTENT(out) :: minv, maxv, mean, mean2, stdev
2437
2438	!!!!!!! Variables
2439	! Nvals: number of values
2440	! vals: values
2441	! minv: minimum value of values
2442	! maxv: maximum value of values
2443	! mean: mean value of values
2444	! mean2: quadratic mean value of values
2445	! stdev: standard deviation of values
2446
2447	minv = MINVAL(vals)
2448	maxv = MAXVAL(vals)
2449
2450	mean=SUM(vals)
2451	mean2=SUM(vals*vals)
2452
2453	mean=mean/Nvals
2454	mean2=mean2/Nvals
2455
2456	stdev=SQRT(mean2 - mean*mean)
2457
2458	RETURN
2459
2460	END SUBROUTINE StatsR_k
2461
2462	END MODULE module_scientific

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

Download in other formats:

Original Format