Context Navigation

module_scientific.f90 @ 1614

Last change on this file since 1614 was 1614, checked in by lfita, 8 years ago
Fixing last issues on polygons detection Adding direct acces to the 'module_scientific.f90' as `module_ForSci'
File size: 31.0 KB

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

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

Download in other formats:

Original Format