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

Rev	Line
[1609]	1	MODULE module_scientific
	2	! Module of the scientific function/subroutines
	3
	4	!!!!!!! Functions & subroutines
[1618]	5	! all_polygons_properties: Subroutine to determine the properties of all polygons in a 2D field:
[1612]	6	! borders_matrixL: Subroutine to provide the borders of a logical array (interested in .TRUE.)
[1650]	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
[1652]	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
[1650]	11	! coincidence_poly: Subtourine to determine which is the coincident polygon when a boolean polygon is provided
	12	! to a map of integer polygons
[1652]	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
[1616]	15	! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
[1609]	16	! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
[1612]	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
[1618]	23	! polygon_properties: Subroutine to determine the properties of a polygon (as .TRUE. matrix)
[1612]	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!
[1650]	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)
[1652]	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
[1609]	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.
[1618]	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
[1609]	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
[1652]	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
[1654]	66	REAL(r_k), DIMENSION(Nmaxpoly) :: Aprevpolys, Acurrpolys
	67	REAL(r_k), DIMENSION(2,Nmaxpoly) :: Cprevpolys, Ccurrpolys
[1652]	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
[1654]	128	Acurrpolys(ip) = areapolys(i,it)
	129	Ccurrpolys(1,ip) = ctrpolys(1,i,it)
	130	Ccurrpolys(2,ip) = ctrpolys(2,i,it)
[1652]	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
[1654]	164	Aprevpolys = Acurrpolys
	165	Cprevpolys = Ccurrpolys
[1652]	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
[1654]	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)
[1652]	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
[1654]	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, &
[1652]	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
[1654]	231	! Which index corresponds to this coincidence?
	232	iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
	233	IF (iindep == -1) THEN
[1652]	234	Nindep = Nindep + 1
	235	SpolysIndep(Nindep,it+1) = coins(i)
	236	END IF
[1654]	237	iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i))
[1652]	238	IF (iindep < 0) THEN
[1654]	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)
[1652]	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.
[1655]	387	finaltracks(2,itrack,it) = SUM(tracks(3,:,itt,it))/Nprev*1.
	388	finaltracks(3,itrack,it) = SUM(tracks(4,:,itt,it))/Nprev*1.
[1652]	389	finaltracks(4,itrack,it) = it*1.
	390	END DO
	391	END DO
	392
[1654]	393	DEALLOCATE(coins)
	394	DEALLOCATE(coinsNpts)
	395
[1652]	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)
[1650]	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
[1651]	579	INTEGER, DIMENSION(dx,dy,dt), INTENT(in) :: allpolys
[1652]	580	REAL(r_k), INTENT(in) :: minarea
[1651]	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), &
[1650]	584	INTENT(out) :: tracks
[1651]	585	REAL(r_k), DIMENSION(4,Nmaxtracks,dt), INTENT(out) :: finaltracks
[1650]	586
	587	! Local
	588	INTEGER :: i, j, ip, it, iip, itt
	589	INTEGER :: ierr
[1651]	590	INTEGER, DIMENSION(Nmaxpoly,dt) :: coincidencies, NOcoincidencies
[1650]	591	INTEGER, DIMENSION(:), ALLOCATABLE :: coins
	592	INTEGER, DIMENSION(:,:), ALLOCATABLE :: coinsNpts
[1651]	593	INTEGER, DIMENSION(Nmaxpoly,dt) :: polycoincidencies
	594	INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt) :: coincidenciesNpts
[1650]	595	INTEGER :: Nmeet, Nsearch, Nindep
	596	INTEGER, DIMENSION(dt) :: Nindeppolys
	597	CHARACTER(len=5) :: NcoinS
[1651]	598	INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt) :: polysIndep
	599	INTEGER, DIMENSION(Nmaxpoly,dt) :: NpolysIndep
	600	INTEGER, DIMENSION(Nmaxpoly,dt) :: SpolysIndep
[1650]	601	INTEGER :: iindep, iiprev
[1651]	602	INTEGER :: Nprev, NNprev, Ntprev
	603	LOGICAL :: Indeppolychained
[1650]	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
[1652]	613	! minarea: minimal area (in same units as areapolys) to perform the tracking
[1650]	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
[1651]	645	SpolysIndep(i,it) = i
	646	NpolysIndep(1:Nmeet,it) = 1
	647	polysIndep(1,i,it) = i
[1650]	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
[1651]	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)
[1650]	667
[1651]	668	polycoincidencies(1:Nmeet,it+1) = coins
	669	coincidenciesNpts(1:Nmeet,1:Nsearch,it+1) = coinsNpts
[1650]	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
[1651]	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
[1650]	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
[1651]	695	IF (coincidencies(ip,it+1) == 0) THEN
[1650]	696	Nindep = Nindep + 1
[1651]	697	SpolysIndep(Nindep,it+1) = ip
[1650]	698	END IF
[1651]	699	coincidencies(ip,it+1) = coincidencies(ip,it+1) + 1
[1650]	700	DO iindep=1, Nindep
[1651]	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
[1650]	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)
[1651]	717	PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', &
	718	polysIndep(i,1:NpolysIndep(i,it+1),it+1)
[1650]	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)
[1651]	729	PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':', &
	730	polysIndep(ip,1:NpolysIndep(ip,it),it)
[1650]	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
[1651]	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
[1650]	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)
[1651]	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)
[1650]	761	Indeppolychained = .FALSE.
	762
	763	! Number of tracks previous time-step
	764	! Looping overall
	765	it1_tracks: DO itt=1, Ntracks(it-1)
[1651]	766	itrack = tracks(1,1,itt,it-1)
[1650]	767	! Number polygons ID assigned
[1651]	768	Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0)
	769	IF (dbg) PRINT *,itt,' track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1)
[1650]	770
	771	! Looking for coincidencies
	772	DO iip=1, Ntprev
[1651]	773	IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN
[1650]	774	Indeppolychained = .TRUE.
[1651]	775	IF (dbg) PRINT *,' coincidence found by polygon:', tracks(2,iip,itt,it-1)
[1650]	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
[1651]	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
[1650]	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
[1651]	802	maxtrack = INT(MAXVAL(tracks(1,:,:,:)*1.))
[1650]	803	IF (dbg) PRINT *,' New track!', maxtrack+1
	804
	805	! Assigning all the IDs to the next step of the track
[1651]	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
[1650]	815	END DO
	816	END IF
	817
	818	END DO current_poly
	819
[1652]	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
[1650]	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.
[1651]	834	maxtrack = MAXVAL(tracks(1,:,:,:))
[1650]	835
	836	DO it=1, dt
	837	DO itt=1, Ntracks(it)
[1651]	838	itrack = INT(tracks(1,1,itt,it))
	839	Nprev = COUNT(INT(tracks(2,:,itt,it)) /= 0)
[1650]	840	PRINT *,'it:', it,'itrack:', itrack, 'Nprev:', Nprev
[1651]	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)
[1650]	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
[1651]	863	REAL(r_k), DIMENSION(2,Nallpoly), INTENT(in) :: icpolys
	864	REAL(r_k), DIMENSION(2,Npoly), INTENT(in) :: cpolys
[1650]	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
[1651]	907	dist = SQRT( (icpolys(1,ip)1.-cpolys(1,j)1.)*2 + (icpolys(2,ip)1.-cpolys(2,j)1.)*2 )
[1650]	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
[1620]	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)
[1618]	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
[1620]	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
[1618]	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
[1620]	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
[1618]	1023
	1024	IMPLICIT NONE
	1025
[1620]	1026	LOGICAL, INTENT(in) :: dbg
[1618]	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
[1620]	1031	CHARACTER(len=1000), INTENT(in) :: projN
[1618]	1032	INTEGER, DIMENSION(Npoly), INTENT(out) :: Npolyptss
[1620]	1033	INTEGER, DIMENSION(Npoly,2), INTENT(out) :: xxtrms, yxtrms, meanctrs
[1618]	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
[1620]	1037	INTEGER, DIMENSION(Npoly,2), INTENT(out) :: nvcoords, xvcoords
	1038	REAL(r_k), DIMENSION(Npoly,2), INTENT(out) :: meanwctrs, meanvnctrs, meanvxctrs
[1618]	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
[1655]	1052	! 'ctsarea': Constant Area
[1618]	1053	! 'lon/lat': for regular longitude-latitude
[1620]	1054	! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
[1618]	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
[1620]	1058	! meanwctrs: lon, lat coordinates of the center from the spatial-weighted mean of the polygons
[1618]	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
[1620]	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
[1618]	1068
	1069	fname = 'all_polygons_properties'
	1070
	1071	! Initializing matrices
	1072	Npolyptss = -1
	1073	xxtrms = fillval64
	1074	yxtrms = fillval64
	1075	meanctrs = fillval64
[1620]	1076	meanwctrs = fillval64
[1618]	1077	areas = fillval64
[1622]	1078	nvals = fillvalI
[1618]	1079	xvals = fillval64
	1080	mvals = fillval64
	1081	m2vals = fillval64
	1082	stdvals = fillval64
	1083	quants = fillval64
[1622]	1084	nvcoords = fillvalI
	1085	xvcoords = fillvalI
[1620]	1086	meanvnctrs = fillval64
	1087	meanvxctrs = fillval64
[1618]	1088
	1089	DO ip=1, Npoly
	1090	boolpoly = polys == ip
[1620]	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,:))
[1618]	1095	END DO
	1096
	1097	RETURN
	1098
	1099	END SUBROUTINE all_polygons_properties
	1100
[1620]	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)
[1618]	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
[1620]	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
[1618]	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
[1620]	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
[1618]	1118
	1119	IMPLICIT NONE
	1120
[1620]	1121	LOGICAL, INTENT(in) :: dbg
[1618]	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
[1620]	1126	CHARACTER(len=1000), INTENT(in) :: projN
[1618]	1127	INTEGER, INTENT(out) :: Npolypts
	1128	INTEGER, DIMENSION(2), INTENT(out) :: xxtrm, yxtrm, meanctr
[1620]	1129	INTEGER, DIMENSION(2), INTENT(out) :: nvcoord, xvcoord
	1130	REAL(r_k), DIMENSION(2), INTENT(out) :: meanwctr, meanvnctr, meanvxctr
[1618]	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
[1620]	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
[1618]	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
[1620]	1155	! 'ctsarea': Constant Area
[1618]	1156	! 'lon/lat': for regular longitude-latitude
[1620]	1157	! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
[1618]	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
[1620]	1161	! meanwctr: lon, lat coordinates of the center from the spatial-weighted mean of the polygon
[1618]	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
[1620]	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
[1618]	1171
	1172	fname = 'polygon_properties'
	1173
[1620]	1174	IF (dbg) PRINT *," '" // TRIM(fname) // "' ..."
	1175
[1618]	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
[1620]	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.
[1654]	1222	xcorr = 0.
	1223	ycorr = 0.
[1620]	1224
	1225	nval = fillval64
	1226	xval = -fillval64
	1227
[1618]	1228	ip = 1
	1229	DO i=1,dx
	1230	DO j=1,dy
	1231	IF (poly(i,j)) THEN
[1622]	1232	polypts(ip,1) = i
	1233	polypts(ip,2) = j
[1618]	1234	polyvals(ip) = values(i,j)
[1620]	1235	SELECT CASE (TRIM(projS))
	1236	CASE ('ctsarea')
	1237	! Constant Area
	1238	xcorr(i,j) = 1.
	1239	ycorr(i,j) = 1.
[1618]	1240	CASE ('lon/lat')
	1241	! Area as fixed yres and sinus-lat varying for xres
[1654]	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
[1620]	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
[1654]	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
[1620]	1256	ycorr(i,j) = 1.
[1618]	1257	END SELECT
[1620]	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
[1618]	1271	ip = ip + 1
	1272	END IF
	1273	END DO
	1274	END DO
	1275
[1620]	1276	IF (dbg) THEN
[1654]	1277	PRINT *,' grid_coord lon lat value _______ '
[1622]	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
[1620]	1282	END IF
	1283
	1284	sumRESx = xres*SUM(xcorr)
	1285	sumRESy = yres*SUM(ycorr)
	1286
[1618]	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 /)
[1620]	1290	meanwctr = (/ meanwctr(1)/sumRESx, meanwctr(2)/sumRESy /)
[1618]	1291
[1620]	1292	IF (dbg) THEN
	1293	PRINT *,' mean grid center: ', meanctr, ' weighted mean center: ', meanwctr
	1294	END IF
[1654]	1295
[1618]	1296	! Statistics of the values within the polygon
	1297	CALL StatsR_K(Npolypts, polyvals, nval, xval, mval, m2val, stdval)
	1298
[1620]	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
[1654]	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
[1620]	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
[1622]	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)
[1620]	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
[1618]	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
[1612]	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
[1614]	1388	Npoly = 0
[1612]	1389
	1390	DO it=1,dt
[1614]	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
[1612]	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
[1614]	1423	LOGICAL, DIMENSION(dx,dy) :: isborder, isbordery
[1612]	1424	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: paths
	1425	INTEGER :: Npath
	1426	INTEGER, DIMENSION(:), ALLOCATABLE :: Nptpaths
	1427	INTEGER, DIMENSION(2) :: xtrx, xtry, meanpth
[1614]	1428	INTEGER :: Nvertx, Npts
[1612]	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
[1618]	1442	! The mathematical maximum woiuld be dx*dy/4, but let's be optimistic... (sorry Jero)
[1614]	1443	Nppt = dx*dy/10
	1444
[1612]	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
[1616]	1460	! Filling with the points of all the space with .TRUE.
	1461	Npts = COUNT(boolmat)
	1462
[1612]	1463	IF (ALLOCATED(points)) DEALLOCATE(points)
[1614]	1464	ALLOCATE(points(Npts,2), STAT=ierr)
[1612]	1465	msg = "Problems allocating matrix 'points'"
	1466	CALL ErrMsg(msg, fname, ierr)
	1467
[1616]	1468	! We only want to localize that points 'inside'
[1612]	1469	ip = 1
	1470	DO i=1, dx
	1471	DO j=1, dy
[1616]	1472	IF (boolmat(i,j)) THEN
	1473	points(ip,1) = i
	1474	points(ip,2) = j
	1475	ip = ip + 1
	1476	END IF
[1612]	1477	END DO
	1478	END DO
	1479
[1614]	1480	CALL borders_matrixL(dx, dy, Nppt, boolmat, borders, isborder, isbordery)
[1612]	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)
[1614]	1492	ALLOCATE(isin(Npts), STAT=ierr)
[1612]	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
[1614]	1514	CALL gridpoints_InsidePolygon(dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), Nvertx,&
	1515	xtrx, xtry, vertxs, Npts, points, isin)
[1612]	1516
	1517	! Filling polygons
[1616]	1518	DO ipp=1, Npts
	1519	IF (isin(ipp)) polys(points(ipp,1),points(ipp,2)) = ip
[1612]	1520	END DO
	1521
[1614]	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
[1612]	1531	END DO
	1532
[1616]	1533	! Cleaning polygons matrix of not-used and paths around holes
	1534	CALL clean_polygons(dx, dy, boolmat, polys, Npoly, dbg)
	1535
[1612]	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
[1616]	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
[1612]	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
[1614]	1768	INTEGER :: Nintersecs, isvertex, ispath
[1612]	1769	INTEGER :: ierr
	1770	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: halo_brdr
[1614]	1771	INTEGER :: Nbrbrdr
[1612]	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
[1614]	1779	! [x/y]pathxtrm extremes of the path
[1612]	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?
[1614]	1809	ispath = index_list_coordsI(Npath, path, (/ix,iy/))
	1810	IF (isbrdr(ix,iy) .AND. (ispath /= -1)) THEN
[1612]	1811	inside(ip) = .TRUE.
	1812	CYCLE
	1813	END IF
	1814
	1815	! Looking along y-axis
[1614]	1816	! Accounting for consecutives borders
	1817	Nbrbrdr = 0
	1818	DO j=MAX(1,ypathxtrm(1)-1),iy-1
[1612]	1819	! Only counting that borders that are not vertexs
[1614]	1820	ispath = index_list_coordsI(Npath, path, (/ix,j/))
[1612]	1821	isvertex = index_list_coordsI(Npath, vrtxs, (/ix,j/))
[1614]	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
[1612]	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
[1614]	1876	lclock = (/ 'W ', 'NW', 'N ', 'NE', 'E ', 'SE', 'S ', 'SW' /)
[1612]	1877	! Transformation to apply
[1614]	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 /)
[1612]	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
[1614]	1912	SUBROUTINE borders_matrixL(dx,dy,dxy,Lmat,brdrs,isbrdr,isbrdry)
[1612]	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
[1614]	1920	LOGICAL, DIMENSION(dx,dy), INTENT(out) :: isbrdr, isbrdry
[1612]	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
[1614]	1931	! isbrdry: matrix with .T./.F. wether the given matrix point is a border or not only along y-axis
[1612]	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
[1614]	1971	isbrdry = isbrdr
	1972
[1612]	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.
[1614]	1995	isbrdry(i,j) = .TRUE.
[1612]	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.
[1614]	2001	isbrdry(i,j+1) = .TRUE.
[1612]	2002	ib=ib+1
	2003	END IF
	2004	END IF
	2005	END DO
	2006	END DO
	2007
[1614]	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
[1612]	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
[1614]	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
[1612]	2164	END IF
[1614]	2165	END DO
	2166	IF (.NOT.found) THEN
[1612]	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
[1609]	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
[1618]	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
[1609]	2462	END MODULE module_scientific

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

Download in other formats:

Original Format