Context Navigation

module_scientific.f90 @ 1617

Last change on this file since 1617 was 1616, checked in by lfita, 8 years ago

Adding:

`clean_polygons': Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole

File size: 34.5 KB

Line
1	MODULE module_scientific
2	! Module of the scientific function/subroutines
3
4	!!!!!!! Functions & subroutines
5	! borders_matrixL: Subroutine to provide the borders of a logical array (interested in .TRUE.)
6	! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
7	! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
8	! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon
9	! following ray casting algorithm
10	! look_clockwise_borders: Subroutine to look clock-wise for a next point within a collection of borders
11	! (limits of a region)
12	! paths_border: Subroutine to search the paths of a border field.
13	! path_properties: Subroutine to determine the properties of a path
14	! polygons: Subroutine to search the polygons of a border field. FORTRAN based. 1st = 1!
15	! polygons_t: Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!
16	! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k)
17	! quantilesR_K: Subroutine to provide the quantiles of a given set of values of type real 'r_k'
18	! rand_sample: Subroutine to randomly sample a range of indices
19	! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order.
20	! SwapR_K*: Subroutine swaps the values of its two formal arguments.
21
22	!!! *Functions/Subroutines to sort values adpated. The method used is usually referred to as "selection" method.
23	! from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
24
25	USE module_definitions
26	USE module_generic
27
28	CONTAINS
29
30	SUBROUTINE polygons_t(dbg, dx, dy, dt, boolmatt, polys, Npoly)
31	! Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!
32
33	IMPLICIT NONE
34
35	INTEGER, INTENT(in) :: dx, dy, dt
36	LOGICAL, DIMENSION(dx,dy,dt), INTENT(in) :: boolmatt
37	LOGICAL, INTENT(in) :: dbg
38	INTEGER, DIMENSION(dt), INTENT(out) :: Npoly
39	INTEGER, DIMENSION(dx,dy,dt), INTENT(out) :: polys
40
41	! Local
42	INTEGER :: i,it
43
44	!!!!!!! Variables
45	! dx,dy: spatial dimensions of the space
46	! boolmatt: boolean matrix tolook for the polygons (.TRUE. based)
47	! polys: found polygons
48	! Npoly: number of polygons found
49
50	fname = 'polygons'
51
52	IF (dbg) PRINT *,TRIM(fname)
53
54	polys = -1
55	Npoly = 0
56
57	DO it=1,dt
58	IF (ANY(boolmatt(:,:,it))) THEN
59	IF (dbg) THEN
60	PRINT *,' it:', it, ' num. TRUE:', COUNT(boolmatt(:,:,it)), 'bool _______'
61	DO i=1,dx
62	PRINT *,boolmatt(i,:,it)
63	END DO
64	END IF
65	CALL polygons(dbg, dx, dy, boolmatt(:,:,it), polys(:,:,it), Npoly(it))
66	ELSE
67	IF (dbg) THEN
68	PRINT *,' it:', it, " without '.TRUE.' values skipiing it!!"
69	END IF
70	END IF
71	END DO
72
73	END SUBROUTINE polygons_t
74
75	SUBROUTINE polygons(dbg, dx, dy, boolmat, polys, Npoly)
76	! Subroutine to search the polygons of a boolean field. FORTRAN based. 1st = 1!
77
78	IMPLICIT NONE
79
80	INTEGER, INTENT(in) :: dx, dy
81	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: boolmat
82	LOGICAL, INTENT(in) :: dbg
83	INTEGER, INTENT(out) :: Npoly
84	INTEGER, DIMENSION(dx,dy), INTENT(out) :: polys
85
86	! Local
87	INTEGER :: i, j, ip, ipp, Nppt
88	INTEGER :: ierr
89	INTEGER, DIMENSION(:,:), ALLOCATABLE :: borders
90	LOGICAL, DIMENSION(dx,dy) :: isborder, isbordery
91	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: paths
92	INTEGER :: Npath
93	INTEGER, DIMENSION(:), ALLOCATABLE :: Nptpaths
94	INTEGER, DIMENSION(2) :: xtrx, xtry, meanpth
95	INTEGER :: Nvertx, Npts
96	INTEGER, DIMENSION(:,:), ALLOCATABLE :: vertxs, points
97	LOGICAL, DIMENSION(:), ALLOCATABLE :: isin
98
99	!!!!!!! Variables
100	! dx,dy: spatial dimensions of the space
101	! boolmat: boolean matrix tolook for the polygons (.TRUE. based)
102	! polys: found polygons
103	! Npoly: number of polygons found
104
105	fname = 'polygons'
106
107	polys = -1
108
109	Nppt = dx*dy/10
110
111	IF (ALLOCATED(borders)) DEALLOCATE(borders)
112	ALLOCATE(borders(Nppt,2), STAT=ierr)
113	msg = "Problems allocating matrix 'borders'"
114	CALL ErrMsg(msg, fname, ierr)
115
116	IF (ALLOCATED(paths)) DEALLOCATE(paths)
117	ALLOCATE(paths(Nppt,Nppt,2), STAT=ierr)
118	msg = "Problems allocating matrix 'paths'"
119	CALL ErrMsg(msg, fname, ierr)
120
121	IF (ALLOCATED(Nptpaths)) DEALLOCATE(Nptpaths)
122	ALLOCATE(Nptpaths(Nppt), STAT=ierr)
123	msg = "Problems allocating matrix 'Nptpaths'"
124	CALL ErrMsg(msg, fname, ierr)
125
126	! Filling with the points of all the space with .TRUE.
127	Npts = COUNT(boolmat)
128
129	IF (ALLOCATED(points)) DEALLOCATE(points)
130	ALLOCATE(points(Npts,2), STAT=ierr)
131	msg = "Problems allocating matrix 'points'"
132	CALL ErrMsg(msg, fname, ierr)
133
134	! We only want to localize that points 'inside'
135	ip = 1
136	DO i=1, dx
137	DO j=1, dy
138	IF (boolmat(i,j)) THEN
139	points(ip,1) = i
140	points(ip,2) = j
141	ip = ip + 1
142	END IF
143	END DO
144	END DO
145
146	CALL borders_matrixL(dx, dy, Nppt, boolmat, borders, isborder, isbordery)
147	CALL paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)
148
149	Npoly = Npath
150
151	DO ip=1, Npath
152	IF (ALLOCATED(vertxs)) DEALLOCATE(vertxs)
153	ALLOCATE(vertxs(Nptpaths(ip),2))
154	msg = "Problems allocating matrix 'vertxs'"
155	CALL ErrMsg(msg, fname, ierr)
156
157	IF (ALLOCATED(isin)) DEALLOCATE(isin)
158	ALLOCATE(isin(Npts), STAT=ierr)
159	msg = "Problems allocating matrix 'isin'"
160	CALL ErrMsg(msg, fname, ierr)
161
162	isin = .FALSE.
163
164	IF (dbg) PRINT *, ' path:', ip, ' N pts:', Nptpaths(ip)
165
166	CALL path_properties(dx, dy, boolmat, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), xtrx, xtry, &
167	meanpth, 'y', Nvertx, vertxs)
168
169	IF (dbg) THEN
170	PRINT *, ' properties _______'
171	PRINT *, ' x-extremes:', xtrx
172	PRINT *, ' y-extremes:', xtry
173	PRINT *, ' center mean:', meanpth
174	PRINT *, ' y-vertexs:', Nvertx,' ________'
175	DO i=1, Nvertx
176	PRINT *,' ',i,':',vertxs(i,:)
177	END DO
178	END IF
179
180	CALL gridpoints_InsidePolygon(dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), Nvertx,&
181	xtrx, xtry, vertxs, Npts, points, isin)
182
183	! Filling polygons
184	DO ipp=1, Npts
185	IF (isin(ipp)) polys(points(ipp,1),points(ipp,2)) = ip
186	END DO
187
188	IF (dbg) THEN
189	PRINT *,' boolmat isborder isbordery polygon (',xtrx(1),',',xtry(1),')x(',xtrx(2),',',xtry(2), &
190	') _______'
191	DO i=xtrx(1), xtrx(2)
192	PRINT *,i,':',boolmat(i,xtry(1):xtry(2)), ' border ', isborder(i,xtry(1):xtry(2)), &
193	' isbordery ', isbordery(i,xtry(1):xtry(2)), ' polygon ', polys(i,xtry(1):xtry(2))
194	END DO
195	END IF
196
197	END DO
198
199	! Cleaning polygons matrix of not-used and paths around holes
200	CALL clean_polygons(dx, dy, boolmat, polys, Npoly, dbg)
201
202	DEALLOCATE (borders)
203	DEALLOCATE (Nptpaths)
204	DEALLOCATE (paths)
205	DEALLOCATE (vertxs)
206	DEALLOCATE (points)
207	DEALLOCATE (isin)
208
209	RETURN
210
211	END SUBROUTINE polygons
212
213	SUBROUTINE clean_polygons(dx, dy, Lmat, pols, Npols, dbg)
214	! Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
215
216	IMPLICIT NONE
217
218	INTEGER, INTENT(in) :: dx, dy
219	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat
220	INTEGER, INTENT(inout) :: Npols
221	INTEGER, DIMENSION(dx,dy), INTENT(inout) :: pols
222	LOGICAL, INTENT(in) :: dbg
223
224	! Local
225	INTEGER :: i,j,ip,iprm
226	INTEGER, DIMENSION(Npols) :: origPol, NotPol, neigPol
227	INTEGER :: ispol, NnotPol
228	CHARACTER(len=4) :: ISa
229
230	!!!!!!! Variables
231	! dx, dy: size of the space
232	! Lmat: original bolean matrix from which the polygons come from
233	! Npols: original number of polygons
234	! pols: polygons space
235
236	fname = 'clean_polygons'
237	IF (dbg) PRINT *," At '" // TRIM(fname) // "' ..."
238
239	origPol = -1
240
241	! Looking for polygons already in space
242	NnotPol = 0
243	DO ip=1, Npols
244	ispol = COUNT(pols-ip == 0)
245	IF (ispol > 0) THEN
246	origPol(ip) = ip
247	ELSE
248	NnotPol = NnotPol + 1
249	NotPol(NnotPol) = ip
250	neigPol(NnotPol) = -1
251	END IF
252	END DO
253
254	IF (dbg) THEN
255	PRINT *,' It should be:', Npols, ' polygons, but already there are:', Npols - NnotPol
256	PRINT *,' Polygons to remove:', NotPol(1:NnotPol)
257	END IF
258
259	! Looking for the hole border of a polygon. This is identify as such polygon point which along
260	! y-axis has NpolygonA, Npolygon, .FALSE.
261	DO i=1,dx
262	DO j=2,dy-1
263	IF ( (pols(i,j-1) /= pols(i,j) .AND. pols(i,j+1) == -1) .AND. (COUNT(NotPol-pols(i,j)==0)==0) &
264	.AND. (pols(i,j) /= -1) .AND. (pols(i,j-1) /= -1)) THEN
265	IF (dbg) PRINT *,' Polygon:', pols(i,j), ' to be removed at point (',i,',',j,'); j-1:', &
266	pols(i,j-1), ' j:', pols(i,j), ' j+1:', pols(i,j+1)
267	NnotPol = NnotPol + 1
268	NotPol(NnotPol) = pols(i,j)
269	neigPol(NnotPol) = pols(i,j-1)
270	END IF
271	END DO
272	END DO
273
274	IF (dbg) THEN
275	PRINT *,' It should be:', Npols, ' polygons, but already there are:', Npols - NnotPol
276	PRINT *,' Polygons to remove after looking for fake border-of-hole polygons _______'
277	DO i=1, NnotPol
278	PRINT *, ' Polygon:', NotPol(i), ' to be replaced by:', neigPol(i)
279	END DO
280	END IF
281
282	! Removing polygons
283	DO iprm=1, NnotPol
284	IF (neigPol(iprm) == -1) THEN
285	WHERE (pols == NotPol(iprm))
286	pols = -1
287	END WHERE
288	IF (dbg) THEN
289	PRINT *,' removing polygon:', NotPol(iprm)
290	END IF
291	ELSE
292	WHERE (pols == NotPol(iprm))
293	pols = neigPol(iprm)
294	END WHERE
295	IF (dbg) THEN
296	PRINT *,' replacing polygon:', NotPol(iprm), ' by:', neigPol(iprm)
297	END IF
298	END IF
299	END DO
300
301	! Re-numbering (descending values)
302	DO i = 1, NnotPol
303	iprm = MAXVAL(NotPol(1:NnotPol))
304	WHERE(pols > iprm)
305	pols = pols - 1
306	END WHERE
307	j = Index1DArrayI(NotPol, NnotPol, iprm)
308	NotPol(j) = -9
309	END DO
310
311	Npols = Npols - NnotPol
312
313	RETURN
314
315	END SUBROUTINE clean_polygons
316
317	SUBROUTINE path_properties(dx, dy, Lmat, Nptspth, pth, xxtrm, yxtrm, meanctr, axs, Nvrtx, vrtxs)
318	! Subroutine to determine the properties of a path:
319	! extremes: minimum and maximum of the path along x,y axes
320	! meancenter: center from the mean of the coordinates of the paths locations
321	! vertexs: path point, without neighbours along a given axis
322
323	IMPLICIT NONE
324
325	INTEGER, INTENT(in) :: dx, dy, Nptspth
326	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat
327	INTEGER, DIMENSION(Nptspth,2), INTENT(in) :: pth
328	CHARACTER, INTENT(in) :: axs
329	INTEGER, DIMENSION(2), INTENT(out) :: meanctr, xxtrm, yxtrm
330	INTEGER, INTENT(out) :: Nvrtx
331	INTEGER, DIMENSION(Nptspth,2), INTENT(out) :: vrtxs
332
333	! Local
334	INTEGER :: i, ip, jp
335	INTEGER :: neig1, neig2
336
337	!!!!!!! Variables
338	! dx,dy: size of the space
339	! Lmat: original matrix of logical values for the path
340	! Nptspth: number of points of the path
341	! pth: path coordinates (clockwise)
342	! axs: axis of finding the vertex
343	! [x/y]xtrm: minimum and maximum coordinates of the path
344	! meanctr: center from the mean of the coordinates of the path
345	! Nvrtx: Number of vertexs of the path
346	! vrtxs: coordinates of the vertexs
347
348	fname = 'path_properties'
349
350	vrtxs = -1
351	Nvrtx = 0
352
353	xxtrm = (/ MINVAL(pth(:,1)), MAXVAL(pth(:,1)) /)
354	yxtrm = (/ MINVAL(pth(:,2)), MAXVAL(pth(:,2)) /)
355	meanctr = (/ SUM(pth(:,1))/Nptspth, SUM(pth(:,2))/Nptspth /)
356
357	IF (axs == 'x' .OR. axs == 'X') THEN
358	! Looking vertexs along x-axis
359	DO i=1, Nptspth
360	ip = pth(i,1)
361	jp = pth(i,2)
362	neig1 = 0
363	neig2 = 0
364	! W-point
365	IF (ip == 1) THEN
366	neig1 = -1
367	ELSE
368	IF (.NOT.Lmat(ip-1,jp)) neig1 = -1
369	END IF
370	! E-point
371	IF (ip == dx) THEN
372	neig2 = -1
373	ELSE
374	IF (.NOT.Lmat(ip+1,jp)) neig2 = -1
375	END IF
376
377	IF (neig1 == -1 .AND. neig2 == -1) THEN
378	Nvrtx = Nvrtx + 1
379	vrtxs(Nvrtx,:) = (/ip,jp/)
380	END IF
381	END DO
382	ELSE IF (axs == 'y' .OR. axs == 'Y') THEN
383	! Looking vertexs along x-axis
384	DO i=1, Nptspth
385	ip = pth(i,1)
386	jp = pth(i,2)
387
388	neig1 = 0
389	neig2 = 0
390	! S-point
391	IF (jp == 1) THEN
392	neig1 = -1
393	ELSE
394	IF (.NOT.Lmat(ip,jp-1)) neig1 = -1
395	END IF
396	! N-point
397	IF (jp == dy) THEN
398	neig2 = -1
399	ELSE
400	IF (.NOT.Lmat(ip,jp+1)) neig2 = -1
401	END IF
402
403	IF (neig1 == -1 .AND. neig2 == -1) THEN
404	Nvrtx = Nvrtx + 1
405	vrtxs(Nvrtx,:) = (/ ip, jp /)
406	END IF
407	END DO
408	ELSE
409	msg = "Axis '" // axs // "' not available" // CHAR(10) // " Available ones: 'x', 'X', 'y, 'Y'"
410	CALL ErrMsg(msg, fname, -1)
411	END IF
412
413	RETURN
414
415	END SUBROUTINE path_properties
416
417	SUBROUTINE gridpoints_InsidePolygon(dx, dy, isbrdr, Npath, path, Nvrtx, xpathxtrm, ypathxtrm, &
418	vrtxs, Npts, pts, inside)
419	! Subroutine to determine if a series of grid points are inside a polygon following ray casting algorithm
420	! FROM: https://en.wikipedia.org/wiki/Point_in_polygon
421
422	IMPLICIT NONE
423
424	INTEGER, INTENT(in) :: dx,dy,Npath,Nvrtx,Npts
425	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isbrdr
426	INTEGER, DIMENSION(Npath,2), INTENT(in) :: path
427	INTEGER, DIMENSION(2), INTENT(in) :: xpathxtrm, ypathxtrm
428	INTEGER, DIMENSION(Npath,2) :: vrtxs
429	INTEGER, DIMENSION(Npts,2), INTENT(in) :: pts
430	LOGICAL, DIMENSION(Npts), INTENT(out) :: inside
431
432	! Local
433	INTEGER :: i,j,ip,ix,iy
434	INTEGER :: Nintersecs, isvertex, ispath
435	INTEGER :: ierr
436	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: halo_brdr
437	INTEGER :: Nbrbrdr
438
439	!!!!!!! Variables
440	! dx,dy: space size
441	! Npath: number of points of the path of the polygon
442	! path: path of the polygon
443	! isbrdr: boolean matrix of the space wqith .T. on polygon border
444	! Nvrtx: number of vertexs of the path
445	! [x/y]pathxtrm extremes of the path
446	! vrtxs: vertexs of the path along y-axis
447	! Npts: number of points
448	! pts: points to look for
449	! inside: vector wether point is inside or not (coincident to a border is inside)
450
451	fname = 'gridpoints_InsidePolygon'
452
453	! Creation of a 1-grid point larger matrix to deal with points reaching the limits
454	IF (ALLOCATED(halo_brdr)) DEALLOCATE(halo_brdr)
455	ALLOCATE(halo_brdr(dx+2,dy+2), STAT=ierr)
456	msg = "Problems allocating matrix 'halo_brdr'"
457	CALL ErrMsg(msg, fname, ierr)
458	halo_brdr = .FALSE.
459
460	DO i=1,dx
461	halo_brdr(i+1,2:dy+1) = isbrdr(i,:)
462	END DO
463
464	inside = .FALSE.
465
466	DO ip=1,Npts
467	Nintersecs = 0
468	ix = pts(ip,1)
469	iy = pts(ip,2)
470	! Point might be outside path range...
471	IF (ix >= xpathxtrm(1) .AND. ix <= xpathxtrm(2) .AND. iy >= ypathxtrm(1) .AND. &
472	iy <= ypathxtrm(2)) THEN
473
474	! It is a border point?
475	ispath = index_list_coordsI(Npath, path, (/ix,iy/))
476	IF (isbrdr(ix,iy) .AND. (ispath /= -1)) THEN
477	inside(ip) = .TRUE.
478	CYCLE
479	END IF
480
481	! Looking along y-axis
482	! Accounting for consecutives borders
483	Nbrbrdr = 0
484	DO j=MAX(1,ypathxtrm(1)-1),iy-1
485	! Only counting that borders that are not vertexs
486	ispath = index_list_coordsI(Npath, path, (/ix,j/))
487	isvertex = index_list_coordsI(Npath, vrtxs, (/ix,j/))
488
489	IF (halo_brdr(ix+1,j+1) .AND. (ispath /= -1) .AND. (isvertex == -1) ) Nintersecs = Nintersecs + 1
490	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
491	Nbrbrdr = Nbrbrdr + 1
492	ELSE
493	! Will remove that consecutive borders above 2
494	IF (Nbrbrdr /= 0) THEN
495	Nintersecs = Nintersecs - MAX(Nbrbrdr-1, 0)
496	Nbrbrdr = 0
497	END IF
498	END IF
499	END DO
500	IF (MOD(Nintersecs,2) /= 0) inside(ip) = .TRUE.
501	END IF
502
503	END DO
504
505	RETURN
506
507	END SUBROUTINE gridpoints_InsidePolygon
508
509	SUBROUTINE look_clockwise_borders(dx,dy,Nbrdrs,brdrs,gbrdr,isbrdr,ix,iy,dbg,xf,yf,iff)
510	! Subroutine to look clock-wise for a next point within a collection of borders (limits of a region)
511
512	IMPLICIT NONE
513
514	INTEGER, INTENT(in) :: dx, dy, Nbrdrs, ix, iy
515	INTEGER, DIMENSION(Nbrdrs,2), INTENT(in) :: brdrs
516	LOGICAL, DIMENSION(Nbrdrs), INTENT(in) :: gbrdr
517	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isbrdr
518	LOGICAL, INTENT(in) :: dbg
519	INTEGER, INTENT(out) :: xf, yf, iff
520
521	! Local
522	INTEGER :: isch
523	CHARACTER(len=2), DIMENSION(8) :: Lclock
524	INTEGER, DIMENSION(8,2) :: spt
525	INTEGER :: iif, jjf
526
527	!!!!!!! Variables
528	! dx, dy: 2D shape ot the space
529	! Nbrdrs: number of brdrs found in this 2D space
530	! brdrs: list of coordinates of the borders
531	! gbrdr: accounts for the use if the given border point
532	! isbrdr: accounts for the matrix of the point is a border or not
533	! ix,iy: coordinates of the point to start to find for
534	! xf,yf: coordinates of the found point
535	! iff: position of the border found within the list of borders
536
537	fname = 'look_clockwise_borders'
538
539	! Looking clock-wise assuming that one starts from the westernmost point
540
541	! Label of the search
542	lclock = (/ 'W ', 'NW', 'N ', 'NE', 'E ', 'SE', 'S ', 'SW' /)
543	! Transformation to apply
544	!spt = (/ (/-1,0/), (/-1,1/), (/0,1/), (/1,1/), (/1,0/), (/1,-1/), (/0,-1/), (/-1,-1/) /)
545	spt(:,1) = (/ -1, -1, 0, 1, 1, 1, 0, -1 /)
546	spt(:,2) = (/ 0, 1, 1, 1, 0, -1, -1, -1 /)
547
548	xf = -1
549	yf = -1
550	DO isch=1, 8
551	! clock-wise search
552	IF (spt(isch,1) >= 0) THEN
553	iif = MIN(dx,ix+spt(isch,1))
554	ELSE
555	iif = MAX(1,ix+spt(isch,1))
556	END IF
557	IF (spt(isch,2) >= 0) THEN
558	jjf = MIN(dy,iy+spt(isch,2))
559	ELSE
560	jjf = MAX(1,iy+spt(isch,2))
561	END IF
562	iff = index_list_coordsI(Nbrdrs, brdrs,(/iif,jjf/))
563	IF (iff > 0) THEN
564	IF (dbg) PRINT *,' ' // lclock(isch) // '-point:', iif,jjf, ':', iff, 'is',isbrdr(iif,jjf), &
565	'got',gbrdr(iff)
566	IF (isbrdr(iif,jjf) .AND. .NOT.gbrdr(iff)) THEN
567	xf = iif
568	yf = jjf
569	EXIT
570	END IF
571	END IF
572	END DO
573
574	RETURN
575
576	END SUBROUTINE look_clockwise_borders
577
578	SUBROUTINE borders_matrixL(dx,dy,dxy,Lmat,brdrs,isbrdr,isbrdry)
579	! Subroutine to provide the borders of a logical array (interested in .TRUE.)
580
581	IMPLICIT NONE
582
583	INTEGER, INTENT(in) :: dx,dy,dxy
584	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat
585	INTEGER, DIMENSION(dxy,2), INTENT(out) :: brdrs
586	LOGICAL, DIMENSION(dx,dy), INTENT(out) :: isbrdr, isbrdry
587
588	! Local
589	INTEGER :: i,j,ib
590
591	!!!!!!! Variables
592	! dx,dy: size of the space
593	! dxy: maximum number of border points
594	! Lmat: Matrix to look for the borders
595	! brdrs: list of coordinates of the borders
596	! isbrdr: matrix with .T./.F. wether the given matrix point is a border or not
597	! isbrdry: matrix with .T./.F. wether the given matrix point is a border or not only along y-axis
598
599	fname = 'borders_matrixL'
600
601	isbrdr = .FALSE.
602	brdrs = -1
603	ib = 1
604
605	! Starting with the borders. If a given point is TRUE it is a path-vertex
606	! Along y-axis
607	DO i=1, dx
608	IF (Lmat(i,1) .AND. .NOT.isbrdr(i,1)) THEN
609	brdrs(ib,1) = i
610	brdrs(ib,2) = 1
611	isbrdr(i,1) = .TRUE.
612	ib=ib+1
613	END IF
614	IF (Lmat(i,dy) .AND. .NOT.isbrdr(i,dy)) THEN
615	brdrs(ib,1) = i
616	brdrs(ib,2) = dy
617	isbrdr(i,dy) = .TRUE.
618	ib=ib+1
619	END IF
620	END DO
621	! Along x-axis
622	DO j=1, dy
623	IF (Lmat(1,j) .AND. .NOT.isbrdr(1,j)) THEN
624	brdrs(ib,1) = 1
625	brdrs(ib,2) = j
626	isbrdr(1,j) = .TRUE.
627	ib=ib+1
628	END IF
629	IF (Lmat(dx,j) .AND. .NOT.isbrdr(dx,j)) THEN
630	brdrs(ib,1) = dx
631	brdrs(ib,2) = j
632	isbrdr(dx,j) = .TRUE.
633	ib=ib+1
634	END IF
635	END DO
636
637	isbrdry = isbrdr
638
639	! Border as that when looking on x-axis points with Lmat(i) /= Lmat(i+1)
640	DO i=1, dx-1
641	DO j=1, dy-1
642	IF ( Lmat(i,j) .NEQV. Lmat(i+1,j) ) THEN
643	IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN
644	brdrs(ib,1) = i
645	brdrs(ib,2) = j
646	isbrdr(i,j) = .TRUE.
647	ib=ib+1
648	ELSE IF (Lmat(i+1,j) .AND. .NOT.isbrdr(i+1,j)) THEN
649	brdrs(ib,1) = i+1
650	brdrs(ib,2) = j
651	isbrdr(i+1,j) = .TRUE.
652	ib=ib+1
653	END IF
654	END IF
655	! y-axis
656	IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN
657	IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN
658	brdrs(ib,1) = i
659	brdrs(ib,2) = j
660	isbrdr(i,j) = .TRUE.
661	isbrdry(i,j) = .TRUE.
662	ib=ib+1
663	ELSE IF (Lmat(i,j+1) .AND. .NOT.isbrdr(i,j+1)) THEN
664	brdrs(ib,1) = i
665	brdrs(ib,2) = j+1
666	isbrdr(i,j+1) = .TRUE.
667	isbrdry(i,j+1) = .TRUE.
668	ib=ib+1
669	END IF
670	END IF
671	END DO
672	END DO
673
674	DO i=1, dx-1
675	DO j=1, dy-1
676	! y-axis
677	IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN
678	IF (Lmat(i,j)) THEN
679	isbrdry(i,j) = .TRUE.
680	ELSE IF (Lmat(i,j+1)) THEN
681	isbrdry(i,j+1) = .TRUE.
682	END IF
683	END IF
684	END DO
685	END DO
686	! only y-axis adding bands of 2 grid points
687	DO i=1, dx-1
688	DO j=2, dy-2
689	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
690	IF (Lmat(i,j)) THEN
691	isbrdry(i,j) = .TRUE.
692	isbrdry(i,j+1) = .TRUE.
693	END IF
694	END IF
695	END DO
696	END DO
697
698	RETURN
699
700	END SUBROUTINE borders_matrixL
701
702	SUBROUTINE paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)
703	! Subroutine to search the paths of a border field.
704
705	IMPLICIT NONE
706
707	INTEGER, INTENT(in) :: dx, dy, Nppt
708	LOGICAL, INTENT(in) :: dbg
709	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isborder
710	INTEGER, DIMENSION(Nppt,2), INTENT(in) :: borders
711	INTEGER, DIMENSION(Nppt,Nppt,2), INTENT(out) :: paths
712	INTEGER, INTENT(out) :: Npath
713	INTEGER, DIMENSION(Nppt), INTENT(out) :: Nptpaths
714
715	! Local
716	INTEGER :: i,j,k,ib
717	INTEGER :: ierr
718	INTEGER :: Nbrdr
719	LOGICAL, DIMENSION(:), ALLOCATABLE :: gotbrdr, emptygotbrdr
720	INTEGER :: iipth, ipath, ip, Nptspath
721	INTEGER :: iib, jjb, iip, ijp, iif, jjf, iff
722	LOGICAL :: found, finishedstep
723
724	!!!!!!! Variables
725	! dx,dy: spatial dimensions of the space
726	! Nppt: possible number of paths and points that the paths can have
727	! isborder: boolean matrix which provide the borders of the polygon
728	! borders: coordinates of the borders of the polygon
729	! paths: coordinates of each found path
730	! Npath: number of paths found
731	! Nptpaths: number of points per path
732
733	fname = 'paths_border'
734
735	! Sarting matrix
736	paths = -1
737	Npath = 0
738	Nptspath = 0
739	Nptpaths = -1
740
741	ib=1
742	finishedstep = .FALSE.
743
744	! Number of border points
745	DO ib=1, Nppt
746	IF (borders(ib,1) == -1 ) EXIT
747	END DO
748	Nbrdr = ib-1
749
750	IF (dbg) THEN
751	PRINT *,' borders _______'
752	DO i=1,Nbrdr
753	PRINT *,' ',i,':',borders(i,:)
754	END DO
755	END IF
756
757	! Matrix which keeps track if a border point has been located
758	IF (ALLOCATED(gotbrdr)) DEALLOCATE(gotbrdr)
759	ALLOCATE(gotbrdr(Nbrdr), STAT=ierr)
760	msg = "Problems allocating matrix 'gotbrdr'"
761	CALL ErrMsg(msg, fname, ierr)
762	IF (ALLOCATED(emptygotbrdr)) DEALLOCATE(emptygotbrdr)
763	ALLOCATE(emptygotbrdr(Nbrdr), STAT=ierr)
764	msg = "Problems allocating matrix 'emptygotbrdr'"
765	CALL ErrMsg(msg, fname, ierr)
766
767	gotbrdr = .FALSE.
768	emptygotbrdr = .FALSE.
769
770	! Starting the fun...
771
772	! Looking along the lines and when a border is found, starting from there in a clock-wise way
773	iipth = 1
774	ipath = 1
775	DO ib=1,Nbrdr
776	iib = borders(iipth,1)
777	jjb = borders(iipth,2)
778	! Starting new path
779	newpath: IF (.NOT.gotbrdr(iipth)) THEN
780	ip = 1
781	Nptspath = 1
782	paths(ipath,ip,:) = borders(iipth,:)
783	gotbrdr(iipth) = .TRUE.
784	! Looking for following clock-wise search
785	! Not looking for W, because search starts from the W
786	iip = iib
787	ijp = jjb
788	DO k=1,Nbrdr
789	IF (dbg) PRINT *,ipath,'iip jip:', iip, ijp
790	found = .FALSE.
791	CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE.,iif,jjf,iff)
792	IF (iif /= -1) THEN
793	ip=ip+1
794	paths(ipath,ip,:) = (/ iif,jjf /)
795	found = .TRUE.
796	gotbrdr(iff) = .TRUE.
797	iip = iif
798	ijp = jjf
799	Nptspath = Nptspath + 1
800	END IF
801
802	IF (dbg) THEN
803	PRINT *,iib,jjb,' end of this round path:', ipath, '_____', gotbrdr
804	DO i=1, Nptspath
805	PRINT *,' ',i,':',paths(ipath,i,:)
806	END DO
807	END IF
808	! If it is not found a next point, might be because it is a non-polygon related value
809	IF (.NOT.found) THEN
810	IF (dbg) PRINT *,'NOT FOUND !!!', gotbrdr
811	! Are still there available borders?
812	IF (ALL(gotbrdr) .EQV. .TRUE.) THEN
813	finishedstep = .TRUE.
814	Npath = ipath
815	Nptpaths(ipath) = Nptspath
816	EXIT
817	ELSE
818	Nptpaths(ipath) = Nptspath
819	! Let's have a look if the previous points in the path have already some 'non-located' neighbourgs
820	DO i=Nptspath,1,-1
821	iip = paths(ipath,i,1)
822	ijp = paths(ipath,i,2)
823	CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE., iif, &
824	jjf,iff)
825	IF (iif /= -1 .AND. iff /= -1) THEN
826	IF (dbg) PRINT *,' re-take path from point:', iif,',',jjf,' n-path:', iff
827	found = .TRUE.
828	iipth = index_list_coordsI(Nppt, borders, (/iip,ijp/))
829	EXIT
830	END IF
831	END DO
832	IF (.NOT.found) THEN
833	! Looking for the next available border point for the new path
834	DO i=1,Nbrdr
835	IF (.NOT.gotbrdr(i)) THEN
836	iipth = i
837	EXIT
838	END IF
839	END DO
840	IF (dbg) PRINT *,' Looking for next path starting at:', iipth, ' point:', &
841	borders(iipth,:)
842	ipath=ipath+1
843	EXIT
844	END IF
845	END IF
846	ELSE
847	IF (dbg) PRINT *,' looking for next point...'
848	END IF
849	IF (finishedstep) EXIT
850	END DO
851	END IF newpath
852	END DO
853	Npath = ipath
854	Nptpaths(ipath) = Nptspath
855
856	DEALLOCATE (gotbrdr)
857	DEALLOCATE (emptygotbrdr)
858
859	RETURN
860
861	END SUBROUTINE paths_border
862
863	SUBROUTINE rand_sample(Nvals, Nsample, sample)
864	! Subroutine to randomly sample a range of indices
865
866	IMPLICIT NONE
867
868	INTEGER, INTENT(in) :: Nvals, Nsample
869	INTEGER, DIMENSION(Nsample), INTENT(out) :: sample
870
871	! Local
872	INTEGER :: i, ind, jmax
873	REAL, DIMENSION(Nsample) :: randv
874	CHARACTER(len=50) :: fname
875	LOGICAL :: found
876	LOGICAL, DIMENSION(Nvals) :: issampled
877	CHARACTER(len=256) :: msg
878	CHARACTER(len=10) :: IS1, IS2
879
880	!!!!!!! Variables
881	! Nvals: number of values
882	! Nsamples: number of samples
883	! sample: samnple
884	fname = 'rand_sample'
885
886	IF (Nsample > Nvals) THEN
887	WRITE(IS1,'(I10)')Nvals
888	WRITE(IS2,'(I10)')Nsample
889	msg = 'Sampling of ' // TRIM(IS1) // ' is too big for ' // TRIM(IS1) // 'values'
890	CALL ErrMsg(msg, fname, -1)
891	END IF
892
893	! Generation of random numbers always the same series during the whole program!
894	CALL RANDOM_NUMBER(randv)
895
896	! Making sure that we do not repeat any value
897	issampled = .FALSE.
898
899	DO i=1, Nsample
900	! Generation of the index from the random numbers
901	ind = MAX(INT(randv(i)*Nvals), 1)
902
903	IF (.NOT.issampled(ind)) THEN
904	sample(i) = ind
905	issampled(ind) = .TRUE.
906	ELSE
907	! Looking around the given index
908	!PRINT *,' Index :', ind, ' already sampled!', issampled(ind)
909	found = .FALSE.
910	DO jmax=1, Nvals
911	ind = MIN(ind+jmax, Nvals)
912	IF (.NOT.issampled(ind)) THEN
913	sample(i) = ind
914	issampled(ind) = .TRUE.
915	found = .TRUE.
916	EXIT
917	END IF
918	ind = MAX(1, ind-jmax)
919	IF (.NOT.issampled(ind)) THEN
920	sample(i) = ind
921	issampled(ind) = .TRUE.
922	found = .TRUE.
923	EXIT
924	END IF
925	END DO
926	IF (.NOT.found) THEN
927	msg = 'sampling could not be finished due to absence of available value!!'
928	CALL ErrMsg(msg, fname, -1)
929	END IF
930	END IF
931
932	END DO
933
934	RETURN
935
936	END SUBROUTINE rand_sample
937
938	SUBROUTINE PrintQuantilesR_K(Nvals, vals, Nquants, qtvs, bspc)
939	! Subroutine to print the quantiles of values REAL(r_k)
940
941	IMPLICIT NONE
942
943	INTEGER, INTENT(in) :: Nvals, Nquants
944	REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals
945	REAL(r_k), DIMENSION(Nquants), INTENT(in) :: qtvs
946	CHARACTER(len=1000), OPTIONAL :: bspc
947
948	! Local
949	INTEGER :: iq
950	LOGICAL, DIMENSION(Nvals) :: search1, search2, search
951	CHARACTER(len=6) :: RS1
952	CHARACTER(len=50) :: fname
953	CHARACTER(len=1000) :: bspcS
954
955	!!!!!!! Variables
956	! vals: series of values
957	! qtvs: values of the quantiles
958	! bspc: base quantity of spaces
959
960	fname = 'PrintQuantilesR_K'
961
962	IF (PRESENT(bspc)) THEN
963	bspcS = bspc
964	ELSE
965	bspcS = ' '
966	END IF
967
968	DO iq=1, Nquants-1
969
970	WHERE (vals >= qtvs(iq))
971	search1 = .TRUE.
972	ELSEWHERE
973	search1 = .FALSE.
974	END WHERE
975
976	WHERE (vals < qtvs(iq+1))
977	search2 = .TRUE.
978	ELSEWHERE
979	search2 = .FALSE.
980	END WHERE
981
982	WHERE (search1 .AND. search2)
983	search = .TRUE.
984	ELSEWHERE
985	search = .FALSE.
986	END WHERE
987
988	WRITE(RS1, '(F6.2)')(iq)*100./(Nquants-1)
989	PRINT *, TRIM(bspcS) // '[',iq,']', TRIM(RS1) // ' %:', qtvs(iq), 'N:', COUNT(search)
990
991	END DO
992
993	RETURN
994
995	END SUBROUTINE PrintQuantilesR_K
996
997	INTEGER FUNCTION FindMinimumR_K(x, dsize, Startv, Endv)
998	! Function returns the location of the minimum in the section between Start and End.
999
1000	IMPLICIT NONE
1001
1002	INTEGER, INTENT(in) :: dsize
1003	REAL(r_k), DIMENSION(dsize), INTENT(in) :: x
1004	INTEGER, INTENT(in) :: Startv, Endv
1005
1006	! Local
1007	REAL(r_k) :: Minimum
1008	INTEGER :: Location
1009	INTEGER :: i
1010
1011	Minimum = x(Startv) ! assume the first is the min
1012	Location = Startv ! record its position
1013	DO i = Startv+1, Endv ! start with next elements
1014	IF (x(i) < Minimum) THEN ! if x(i) less than the min?
1015	Minimum = x(i) ! Yes, a new minimum found
1016	Location = i ! record its position
1017	END IF
1018	END DO
1019
1020	FindMinimumR_K = Location ! return the position
1021
1022	END FUNCTION FindMinimumR_K
1023
1024	SUBROUTINE SwapR_K(a, b)
1025	! Subroutine swaps the values of its two formal arguments.
1026
1027	IMPLICIT NONE
1028
1029	REAL(r_k), INTENT(INOUT) :: a, b
1030	! Local
1031	REAL(r_k) :: Temp
1032
1033	Temp = a
1034	a = b
1035	b = Temp
1036
1037	END SUBROUTINE SwapR_K
1038
1039	SUBROUTINE SortR_K(x, Nx)
1040	! Subroutine receives an array x() r_K and sorts it into ascending order.
1041
1042	IMPLICIT NONE
1043
1044	INTEGER, INTENT(IN) :: Nx
1045	REAL(r_k), DIMENSION(Nx), INTENT(INOUT) :: x
1046
1047	! Local
1048	INTEGER :: i
1049	INTEGER :: Location
1050
1051	DO i = 1, Nx-1 ! except for the last
1052	Location = FindMinimumR_K(x, Nx-i+1, i, Nx) ! find min from this to last
1053	CALL SwapR_K(x(i), x(Location)) ! swap this and the minimum
1054	END DO
1055
1056	END SUBROUTINE SortR_K
1057
1058	SUBROUTINE quantilesR_K(Nvals, vals, Nquants, quants)
1059	! Subroutine to provide the quantiles of a given set of values of type real 'r_k'
1060
1061	IMPLICIT NONE
1062
1063	INTEGER, INTENT(in) :: Nvals, Nquants
1064	REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals
1065	REAL(r_k), DIMENSION(Nquants), INTENT(out) :: quants
1066
1067	! Local
1068	INTEGER :: i
1069	REAL(r_k) :: minv, maxv
1070	REAL(r_k), DIMENSION(Nvals) :: sortedvals
1071
1072	!!!!!!! Variables
1073	! Nvals: number of values
1074	! Rk: kind of real of the values
1075	! vals: values
1076	! Nquants: number of quants
1077	! quants: values at which the quantile start
1078
1079	minv = MINVAL(vals)
1080	maxv = MAXVAL(vals)
1081
1082	sortedvals = vals
1083	! Using from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
1084	CALL SortR_K(sortedvals, Nvals)
1085
1086	quants(1) = minv
1087	DO i=2, Nquants
1088	quants(i) = sortedvals(INT((i-1)*Nvals/Nquants))
1089	END DO
1090
1091	END SUBROUTINE quantilesR_K
1092
1093	END MODULE module_scientific

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

Download in other formats:

Original Format