Context Navigation

module_scientific.f90 @ 1638

Last change on this file since 1638 was 1622, checked in by lfita, 8 years ago
Fixing various issues
File size: 49.3 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.)
[1616]	7	! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
[1609]	8	! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
[1612]	9	! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon
	10	! following ray casting algorithm
	11	! look_clockwise_borders: Subroutine to look clock-wise for a next point within a collection of borders
	12	! (limits of a region)
	13	! paths_border: Subroutine to search the paths of a border field.
	14	! path_properties: Subroutine to determine the properties of a path
[1618]	15	! polygon_properties: Subroutine to determine the properties of a polygon (as .TRUE. matrix)
[1612]	16	! polygons: Subroutine to search the polygons of a border field. FORTRAN based. 1st = 1!
	17	! polygons_t: Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!
[1609]	18	! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k)
	19	! quantilesR_K: Subroutine to provide the quantiles of a given set of values of type real 'r_k'
	20	! rand_sample: Subroutine to randomly sample a range of indices
	21	! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order.
[1618]	22	! StatsR_K: Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a
	23	! series of r_k numbers
[1609]	24	! SwapR_K*: Subroutine swaps the values of its two formal arguments.
	25
	26	!!! *Functions/Subroutines to sort values adpated. The method used is usually referred to as "selection" method.
	27	! from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
	28
	29	USE module_definitions
	30	USE module_generic
	31
	32	CONTAINS
	33
[1620]	34	SUBROUTINE all_polygons_properties(dbg, dx, dy, Npoly, polys, lon, lat, values, xres, yres, projN, &
	35	Npolyptss, xxtrms, yxtrms, meanctrs, meanwctrs, areas, nvals, xvals, mvals, m2vals, stdvals, &
	36	Nquant, quants, nvcoords, xvcoords, meanvnctrs, meanvxctrs)
[1618]	37	! Subroutine to determine the properties of all polygons in a 2D field:
	38	! Number of grid points
	39	! grid-point coordinates of the minimum and maximum of the path along x,y axes
[1620]	40	! grid coordinates of center from the mean of the coordinates of the poygon locations
	41	! lon, lat center from the area weighted mean of the coordinates of the polygon locations
[1618]	42	! area of the polygon (km2)
	43	! minimum and maximum of the values within the polygon
	44	! mean of the values within the polygon
	45	! quadratic mean of the values within the polygon
	46	! standard deviation of the values within the polygon
	47	! number of quantiles
	48	! quantiles of the values within the polygon
[1620]	49	! grid coordinates of the minimum, maximum value within the polygon
	50	! lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon
[1618]	51
	52	IMPLICIT NONE
	53
[1620]	54	LOGICAL, INTENT(in) :: dbg
[1618]	55	INTEGER, INTENT(in) :: dx, dy, Npoly, Nquant
	56	INTEGER, DIMENSION(dx,dy), INTENT(in) :: polys
	57	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: lon, lat, values
	58	REAL(r_k), INTENT(in) :: xres, yres
[1620]	59	CHARACTER(len=1000), INTENT(in) :: projN
[1618]	60	INTEGER, DIMENSION(Npoly), INTENT(out) :: Npolyptss
[1620]	61	INTEGER, DIMENSION(Npoly,2), INTENT(out) :: xxtrms, yxtrms, meanctrs
[1618]	62	REAL(r_k), DIMENSION(Npoly), INTENT(out) :: areas
	63	REAL(r_k), DIMENSION(Npoly), INTENT(out) :: nvals, xvals, mvals, m2vals, stdvals
	64	REAL(r_k), DIMENSION(Npoly, Nquant), INTENT(out) :: quants
[1620]	65	INTEGER, DIMENSION(Npoly,2), INTENT(out) :: nvcoords, xvcoords
	66	REAL(r_k), DIMENSION(Npoly,2), INTENT(out) :: meanwctrs, meanvnctrs, meanvxctrs
[1618]	67
	68	! Local
	69	INTEGER :: ip
	70	LOGICAL, DIMENSION(dx,dy) :: boolpoly
	71
	72	!!!!!!! Variables
	73	! dx,dy: size of the space
	74	! Npoly: number of polygons
	75	! polys: polygon matrix with all polygons (as integer number per polygon)
	76	! lon, lat: geographical coordinates of the grid-points of the matrix
	77	! values: values of the 2D field to use
	78	! [x/y]res resolution along the x and y axis
	79	! projN: name of the projection
	80	! 'lon/lat': for regular longitude-latitude
[1620]	81	! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
[1618]	82	! Npolyptss: number of points of the polygons
	83	! [x/y]xtrms: grid-point coordinates of minimum and maximum coordinates of the polygons
	84	! meanctrs: center from the mean of the coordinates of the polygons
[1620]	85	! meanwctrs: lon, lat coordinates of the center from the spatial-weighted mean of the polygons
[1618]	86	! areas: area of the polygons [km]
	87	! [n/x]vals: minimum and maximum of the values within the polygons
	88	! mvals: mean of the values within the polygons
	89	! m2vals: quadratic mean of the values within the polygons
	90	! stdvals: standard deviation of the values within the polygons
	91	! Nquant: number of quantiles
	92	! quants: quantiles of the values within the polygons
[1620]	93	! [n/x]vcoords: grid coordinates of the minimum/maximum of the values within the polygons
	94	! 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]	95
	96	fname = 'all_polygons_properties'
	97
	98	! Initializing matrices
	99	Npolyptss = -1
	100	xxtrms = fillval64
	101	yxtrms = fillval64
	102	meanctrs = fillval64
[1620]	103	meanwctrs = fillval64
[1618]	104	areas = fillval64
[1622]	105	nvals = fillvalI
[1618]	106	xvals = fillval64
	107	mvals = fillval64
	108	m2vals = fillval64
	109	stdvals = fillval64
	110	quants = fillval64
[1622]	111	nvcoords = fillvalI
	112	xvcoords = fillvalI
[1620]	113	meanvnctrs = fillval64
	114	meanvxctrs = fillval64
[1618]	115
	116	DO ip=1, Npoly
	117	boolpoly = polys == ip
[1620]	118	CALL polygon_properties(dbg, dx, dy, boolpoly, lon, lat, values, xres, yres, projN, Npolyptss(ip),&
	119	xxtrms(ip,:), yxtrms(ip,:), meanctrs(ip,:), meanwctrs(ip,:), areas(ip), nvals(ip), xvals(ip), &
	120	mvals(ip), m2vals(ip), stdvals(ip), Nquant, quants(ip,:), nvcoords(ip,:), xvcoords(ip,:), &
	121	meanvnctrs(ip,:), meanvxctrs(ip,:))
[1618]	122	END DO
	123
	124	RETURN
	125
	126	END SUBROUTINE all_polygons_properties
	127
[1620]	128	SUBROUTINE polygon_properties(dbg, dx, dy, poly, lon, lat, values, xres, yres, projN, Npolypts, &
	129	xxtrm, yxtrm, meanctr, meanwctr, area, nval, xval, mval, m2val, stdval, Nquant, quant, nvcoord, &
	130	xvcoord, meanvnctr, meanvxctr)
[1618]	131	! Subroutine to determine the properties of a polygon (as .TRUE. matrix)
	132	! Number of grid points
	133	! grid-point coordinates of the minimum and maximum of the path along x,y axes
[1620]	134	! grid coordinates of center from the mean of the coordinates of the poygon locations
	135	! lon, lat center from the area weighted mean of the coordinates of the polygon locations
[1618]	136	! area of the polygon (km2)
	137	! minimum and maximum of the values within the polygon
	138	! mean of the values within the polygon
	139	! quadratic mean of the values within the polygon
	140	! standard deviation of the values within the polygon
	141	! number of quantiles
	142	! quantiles of the values within the polygon
[1620]	143	! grid coordinates of the minimum, maximum value within the polygon
	144	! lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon
[1618]	145
	146	IMPLICIT NONE
	147
[1620]	148	LOGICAL, INTENT(in) :: dbg
[1618]	149	INTEGER, INTENT(in) :: dx, dy, Nquant
	150	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: poly
	151	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: lon, lat, values
	152	REAL(r_k), INTENT(in) :: xres, yres
[1620]	153	CHARACTER(len=1000), INTENT(in) :: projN
[1618]	154	INTEGER, INTENT(out) :: Npolypts
	155	INTEGER, DIMENSION(2), INTENT(out) :: xxtrm, yxtrm, meanctr
[1620]	156	INTEGER, DIMENSION(2), INTENT(out) :: nvcoord, xvcoord
	157	REAL(r_k), DIMENSION(2), INTENT(out) :: meanwctr, meanvnctr, meanvxctr
[1618]	158	REAL(r_k), INTENT(out) :: area
	159	REAL(r_k), INTENT(out) :: nval, xval, mval, m2val, stdval
	160	REAL(r_k), DIMENSION(Nquant), INTENT(out) :: quant
	161
	162	! Local
	163	INTEGER :: i, j, ip
	164	INTEGER :: ierr
	165	INTEGER, DIMENSION(:,:), ALLOCATABLE :: polypts
[1620]	166	REAL(r_k), DIMENSION(:), ALLOCATABLE :: polyvals, distvn, distvx
	167	REAL(r_k) :: lonNADIR, latNADIR
	168	REAL(r_k) :: sumRESx, sumRESy
	169	REAL(r_k), DIMENSION(dx,dy) :: xcorr, ycorr
	170	CHARACTER(len=200), DIMENSION(3) :: satSvals
	171	CHARACTER(len=50) :: projS
	172	REAL(r_k) :: sumDISTnlon, sumDISTnlat, sumDISTxlon, &
	173	sumDISTxlat
[1618]	174
	175	!!!!!!! Variables
	176	! dx,dy: size of the space
	177	! poly: polygon matrix (boolean)
	178	! lon, lat: geographical coordinates of the grid-points of the matrix
	179	! values: values of the 2D field to use
	180	! [x/y]res resolution along the x and y axis
	181	! projN: name of the projection
[1620]	182	! 'ctsarea': Constant Area
[1618]	183	! 'lon/lat': for regular longitude-latitude
[1620]	184	! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
[1618]	185	! Npolypts: number of points of the polygon
	186	! [x/y]xtrm: grid-point coordinates of minimum and maximum coordinates of the polygon
	187	! meanctr: center from the mean of the coordinates of the polygon
[1620]	188	! meanwctr: lon, lat coordinates of the center from the spatial-weighted mean of the polygon
[1618]	189	! area: area of the polygon [km]
	190	! [n/x]val: minimum and maximum of the values within the polygon
	191	! mval: mean of the values within the polygon
	192	! m2val: quadratic mean of the values within the polygon
	193	! stdval: standard deviation of the values within the polygon
	194	! Nquant: number of quantiles
[1620]	195	! quant: quantiles of the values within the polygon
	196	! [n/x]vcoord: grid coordinates of the minimum/maximum of the values within the polygon
	197	! 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]	198
	199	fname = 'polygon_properties'
	200
[1620]	201	IF (dbg) PRINT *," '" // TRIM(fname) // "' ..."
	202
[1618]	203	! Getting grid-point coordinates of the polygon
	204	Npolypts = COUNT(poly)
	205
	206	IF (ALLOCATED(polypts)) DEALLOCATE(polypts)
	207	ALLOCATE(polypts(Npolypts,2), STAT=ierr)
	208	msg = "Problems allocating 'polypts'"
	209	CALL ErrMsg(msg, fname, ierr)
	210
	211	IF (ALLOCATED(polyvals)) DEALLOCATE(polyvals)
	212	ALLOCATE(polyvals(Npolypts), STAT=ierr)
	213	msg = "Problems allocating 'polyvals'"
	214	CALL ErrMsg(msg, fname, ierr)
	215
[1620]	216	IF (ALLOCATED(distvn)) DEALLOCATE(distvn)
	217	ALLOCATE(distvn(Npolypts), STAT=ierr)
	218	msg = "Problems allocating 'distvn'"
	219	CALL ErrMsg(msg, fname, ierr)
	220
	221	IF (ALLOCATED(distvx)) DEALLOCATE(distvx)
	222	ALLOCATE(distvx(Npolypts), STAT=ierr)
	223	msg = "Problems allocating 'distvx'"
	224	CALL ErrMsg(msg, fname, ierr)
	225
	226	IF (projN(1:7) == 'lon/lat') THEN
	227	projS = projN
	228	ELSE IF (projN(1:7) == 'ctsarea') THEN
	229	projS = projN
	230	ELSE IF (projN(1:9) == 'nadir-sat') THEN
	231	projS = 'nadir-sat'
	232	CALL split(projN, ',', 3, satSvals)
	233	READ(satSvals(2),'(F200.100)')lonNadir
	234	READ(satSvals(3),'(F200.100)')latNadir
	235	IF (dbg) PRINT *," 'nadir-geostationary-satellite' based projection of data with nadir " // &
	236	"location at:", lonNadir, latNadir
	237	ELSE
	238	msg = "Projection '" // TRIM(projN) // "' not ready" // CHAR(10) // " available ones: " // &
	239	"'ctsarea', 'lon/lat', 'nadir-sat'"
	240	CALL ErrMsg(msg,fname,-1)
	241	END IF
	242
	243	area = 0.
	244	sumRESx = 0.
	245	sumRESy = 0.
	246	meanwctr = 0.
	247	meanvnctr = 0.
	248	meanvxctr = 0.
	249
	250	nval = fillval64
	251	xval = -fillval64
	252
[1618]	253	ip = 1
	254	DO i=1,dx
	255	DO j=1,dy
	256	IF (poly(i,j)) THEN
[1622]	257	polypts(ip,1) = i
	258	polypts(ip,2) = j
[1618]	259	polyvals(ip) = values(i,j)
[1620]	260	SELECT CASE (TRIM(projS))
	261	CASE ('ctsarea')
	262	! Constant Area
	263	xcorr(i,j) = 1.
	264	ycorr(i,j) = 1.
[1618]	265	CASE ('lon/lat')
	266	! Area as fixed yres and sinus-lat varying for xres
[1620]	267	xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
	268	ycorr(i,j) = 1.
	269	CASE ('nadir-sat')
	270	! Area from nadir resolution and degrading as we get far from satellite's nadir
	271	! GOES-E: 0 N, 75 W
	272	xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
	273	ycorr(i,j) = 1.
[1618]	274	END SELECT
[1620]	275	area = area + xresxcorr(i,j)yres*ycorr(i,j)
	276	meanwctr(1) = meanwctr(1) + lon(i,j)xresxcorr(i,j)
	277	meanwctr(2) = meanwctr(2) + lat(i,j)yresycorr(i,j)
	278	IF (nval > values(i,j)) THEN
	279	nvcoord(1) = i
	280	nvcoord(2) = j
	281	nval = values(i,j)
	282	END IF
	283	IF (xval < values(i,j)) THEN
	284	xvcoord(1) = i
	285	xvcoord(2) = j
	286	xval = values(i,j)
	287	END IF
[1618]	288	ip = ip + 1
	289	END IF
	290	END DO
	291	END DO
	292
[1620]	293	IF (dbg) THEN
[1622]	294	PRINT *,' grdi_coord lon lat value _______ '
	295	DO ip=1, Npolypts
	296	PRINT *, polypts(ip,:), ';', lon(polypts(ip,1),polypts(ip,2)), lat(polypts(ip,1),polypts(ip,2)),&
	297	':', polyvals(ip)
	298	END DO
[1620]	299	END IF
	300
	301	sumRESx = xres*SUM(xcorr)
	302	sumRESy = yres*SUM(ycorr)
	303
[1618]	304	xxtrm = (/ MINVAL(polypts(:,1)), MAXVAL(polypts(:,1)) /)
	305	yxtrm = (/ MINVAL(polypts(:,2)), MAXVAL(polypts(:,2)) /)
	306	meanctr = (/ SUM(polypts(:,1))/Npolypts, SUM(polypts(:,2))/Npolypts /)
[1620]	307	meanwctr = (/ meanwctr(1)/sumRESx, meanwctr(2)/sumRESy /)
[1618]	308
[1620]	309	IF (dbg) THEN
	310	PRINT *,' mean grid center: ', meanctr, ' weighted mean center: ', meanwctr
	311	END IF
	312
[1618]	313	! Statistics of the values within the polygon
	314	CALL StatsR_K(Npolypts, polyvals, nval, xval, mval, m2val, stdval)
	315
[1620]	316	IF (dbg) THEN
	317	PRINT *,' minimum value: ', nval, ' maximum:', xval, ' mean:', mval
	318	PRINT *,' coor. minimum: ', nvcoord
	319	PRINT *,' coor. maximum: ', xvcoord
	320	END IF
	321
	322	! Mean center weighted to minimum and maximum values
	323	distvn = DABS(polyvals - nval)
	324	distvx = DABS(polyvals - xval)
	325
	326	meanvnctr = 0.
	327	meanvxctr = 0.
	328	sumDISTnlon = 0.
	329	sumDISTnlat = 0.
	330	sumDISTxlon = 0.
	331	sumDISTxlat = 0.
	332
	333	ip = 1
	334	DO i=1,dx
	335	DO j=1,dy
	336	IF (poly(i,j)) THEN
	337	meanvnctr(1) = meanvnctr(1)+lon(i,j)distvn(ip)xres*xcorr(i,j)
	338	meanvnctr(2) = meanvnctr(2)+lat(i,j)distvn(ip)yres*ycorr(i,j)
	339
	340	meanvxctr(1) = meanvxctr(1)+lon(i,j)distvx(ip)xres*xcorr(i,j)
	341	meanvxctr(2) = meanvxctr(2)+lat(i,j)distvx(ip)yres*ycorr(i,j)
	342
[1622]	343	sumDISTnlon = sumDISTnlon + distvn(ip)xresxcorr(i,j)
	344	sumDISTnlat = sumDISTnlat + distvn(ip)yresycorr(i,j)
	345	sumDISTxlon = sumDISTxlon + distvx(ip)xresxcorr(i,j)
	346	sumDISTxlat = sumDISTxlat + distvx(ip)yresycorr(i,j)
[1620]	347
	348	ip = ip + 1
	349	END IF
	350	END DO
	351	END DO
	352
	353	meanvnctr = (/ meanvnctr(1)/sumDISTnlon, meanvnctr(2)/sumDISTnlat /)
	354	meanvxctr = (/ meanvxctr(1)/sumDISTxlon, meanvxctr(2)/sumDISTxlat /)
	355
	356	IF (dbg) THEN
	357	PRINT *,' mean center to minimum: ', meanvnctr, ' to maximum: ', meanvxctr
	358	END IF
	359
[1618]	360	! Quantiles of the values within the polygon
	361	quant = -9999.d0
	362	IF (Npolypts > Nquant) THEN
	363	CALL quantilesR_K(Npolypts, polyvals, Nquant, quant)
	364	ELSE
	365	CALL SortR_K(polyvals, Npolypts)
	366	END IF
	367
	368	DEALLOCATE (polypts)
	369	DEALLOCATE (polyvals)
	370
	371	RETURN
	372
	373	END SUBROUTINE polygon_properties
	374
[1612]	375	SUBROUTINE polygons_t(dbg, dx, dy, dt, boolmatt, polys, Npoly)
	376	! Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!
	377
	378	IMPLICIT NONE
	379
	380	INTEGER, INTENT(in) :: dx, dy, dt
	381	LOGICAL, DIMENSION(dx,dy,dt), INTENT(in) :: boolmatt
	382	LOGICAL, INTENT(in) :: dbg
	383	INTEGER, DIMENSION(dt), INTENT(out) :: Npoly
	384	INTEGER, DIMENSION(dx,dy,dt), INTENT(out) :: polys
	385
	386	! Local
	387	INTEGER :: i,it
	388
	389	!!!!!!! Variables
	390	! dx,dy: spatial dimensions of the space
	391	! boolmatt: boolean matrix tolook for the polygons (.TRUE. based)
	392	! polys: found polygons
	393	! Npoly: number of polygons found
	394
	395	fname = 'polygons'
	396
	397	IF (dbg) PRINT *,TRIM(fname)
	398
	399	polys = -1
[1614]	400	Npoly = 0
[1612]	401
	402	DO it=1,dt
[1614]	403	IF (ANY(boolmatt(:,:,it))) THEN
	404	IF (dbg) THEN
	405	PRINT *,' it:', it, ' num. TRUE:', COUNT(boolmatt(:,:,it)), 'bool _______'
	406	DO i=1,dx
	407	PRINT *,boolmatt(i,:,it)
	408	END DO
	409	END IF
	410	CALL polygons(dbg, dx, dy, boolmatt(:,:,it), polys(:,:,it), Npoly(it))
	411	ELSE
	412	IF (dbg) THEN
	413	PRINT *,' it:', it, " without '.TRUE.' values skipiing it!!"
	414	END IF
[1612]	415	END IF
	416	END DO
	417
	418	END SUBROUTINE polygons_t
	419
	420	SUBROUTINE polygons(dbg, dx, dy, boolmat, polys, Npoly)
	421	! Subroutine to search the polygons of a boolean field. FORTRAN based. 1st = 1!
	422
	423	IMPLICIT NONE
	424
	425	INTEGER, INTENT(in) :: dx, dy
	426	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: boolmat
	427	LOGICAL, INTENT(in) :: dbg
	428	INTEGER, INTENT(out) :: Npoly
	429	INTEGER, DIMENSION(dx,dy), INTENT(out) :: polys
	430
	431	! Local
	432	INTEGER :: i, j, ip, ipp, Nppt
	433	INTEGER :: ierr
	434	INTEGER, DIMENSION(:,:), ALLOCATABLE :: borders
[1614]	435	LOGICAL, DIMENSION(dx,dy) :: isborder, isbordery
[1612]	436	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: paths
	437	INTEGER :: Npath
	438	INTEGER, DIMENSION(:), ALLOCATABLE :: Nptpaths
	439	INTEGER, DIMENSION(2) :: xtrx, xtry, meanpth
[1614]	440	INTEGER :: Nvertx, Npts
[1612]	441	INTEGER, DIMENSION(:,:), ALLOCATABLE :: vertxs, points
	442	LOGICAL, DIMENSION(:), ALLOCATABLE :: isin
	443
	444	!!!!!!! Variables
	445	! dx,dy: spatial dimensions of the space
	446	! boolmat: boolean matrix tolook for the polygons (.TRUE. based)
	447	! polys: found polygons
	448	! Npoly: number of polygons found
	449
	450	fname = 'polygons'
	451
	452	polys = -1
	453
[1618]	454	! The mathematical maximum woiuld be dx*dy/4, but let's be optimistic... (sorry Jero)
[1614]	455	Nppt = dx*dy/10
	456
[1612]	457	IF (ALLOCATED(borders)) DEALLOCATE(borders)
	458	ALLOCATE(borders(Nppt,2), STAT=ierr)
	459	msg = "Problems allocating matrix 'borders'"
	460	CALL ErrMsg(msg, fname, ierr)
	461
	462	IF (ALLOCATED(paths)) DEALLOCATE(paths)
	463	ALLOCATE(paths(Nppt,Nppt,2), STAT=ierr)
	464	msg = "Problems allocating matrix 'paths'"
	465	CALL ErrMsg(msg, fname, ierr)
	466
	467	IF (ALLOCATED(Nptpaths)) DEALLOCATE(Nptpaths)
	468	ALLOCATE(Nptpaths(Nppt), STAT=ierr)
	469	msg = "Problems allocating matrix 'Nptpaths'"
	470	CALL ErrMsg(msg, fname, ierr)
	471
[1616]	472	! Filling with the points of all the space with .TRUE.
	473	Npts = COUNT(boolmat)
	474
[1612]	475	IF (ALLOCATED(points)) DEALLOCATE(points)
[1614]	476	ALLOCATE(points(Npts,2), STAT=ierr)
[1612]	477	msg = "Problems allocating matrix 'points'"
	478	CALL ErrMsg(msg, fname, ierr)
	479
[1616]	480	! We only want to localize that points 'inside'
[1612]	481	ip = 1
	482	DO i=1, dx
	483	DO j=1, dy
[1616]	484	IF (boolmat(i,j)) THEN
	485	points(ip,1) = i
	486	points(ip,2) = j
	487	ip = ip + 1
	488	END IF
[1612]	489	END DO
	490	END DO
	491
[1614]	492	CALL borders_matrixL(dx, dy, Nppt, boolmat, borders, isborder, isbordery)
[1612]	493	CALL paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)
	494
	495	Npoly = Npath
	496
	497	DO ip=1, Npath
	498	IF (ALLOCATED(vertxs)) DEALLOCATE(vertxs)
	499	ALLOCATE(vertxs(Nptpaths(ip),2))
	500	msg = "Problems allocating matrix 'vertxs'"
	501	CALL ErrMsg(msg, fname, ierr)
	502
	503	IF (ALLOCATED(isin)) DEALLOCATE(isin)
[1614]	504	ALLOCATE(isin(Npts), STAT=ierr)
[1612]	505	msg = "Problems allocating matrix 'isin'"
	506	CALL ErrMsg(msg, fname, ierr)
	507
	508	isin = .FALSE.
	509
	510	IF (dbg) PRINT *, ' path:', ip, ' N pts:', Nptpaths(ip)
	511
	512	CALL path_properties(dx, dy, boolmat, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), xtrx, xtry, &
	513	meanpth, 'y', Nvertx, vertxs)
	514
	515	IF (dbg) THEN
	516	PRINT *, ' properties _______'
	517	PRINT *, ' x-extremes:', xtrx
	518	PRINT *, ' y-extremes:', xtry
	519	PRINT *, ' center mean:', meanpth
	520	PRINT *, ' y-vertexs:', Nvertx,' ________'
	521	DO i=1, Nvertx
	522	PRINT *,' ',i,':',vertxs(i,:)
	523	END DO
	524	END IF
	525
[1614]	526	CALL gridpoints_InsidePolygon(dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), Nvertx,&
	527	xtrx, xtry, vertxs, Npts, points, isin)
[1612]	528
	529	! Filling polygons
[1616]	530	DO ipp=1, Npts
	531	IF (isin(ipp)) polys(points(ipp,1),points(ipp,2)) = ip
[1612]	532	END DO
	533
[1614]	534	IF (dbg) THEN
	535	PRINT *,' boolmat isborder isbordery polygon (',xtrx(1),',',xtry(1),')x(',xtrx(2),',',xtry(2), &
	536	') _______'
	537	DO i=xtrx(1), xtrx(2)
	538	PRINT *,i,':',boolmat(i,xtry(1):xtry(2)), ' border ', isborder(i,xtry(1):xtry(2)), &
	539	' isbordery ', isbordery(i,xtry(1):xtry(2)), ' polygon ', polys(i,xtry(1):xtry(2))
	540	END DO
	541	END IF
	542
[1612]	543	END DO
	544
[1616]	545	! Cleaning polygons matrix of not-used and paths around holes
	546	CALL clean_polygons(dx, dy, boolmat, polys, Npoly, dbg)
	547
[1612]	548	DEALLOCATE (borders)
	549	DEALLOCATE (Nptpaths)
	550	DEALLOCATE (paths)
	551	DEALLOCATE (vertxs)
	552	DEALLOCATE (points)
	553	DEALLOCATE (isin)
	554
	555	RETURN
	556
	557	END SUBROUTINE polygons
	558
[1616]	559	SUBROUTINE clean_polygons(dx, dy, Lmat, pols, Npols, dbg)
	560	! Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
	561
	562	IMPLICIT NONE
	563
	564	INTEGER, INTENT(in) :: dx, dy
	565	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat
	566	INTEGER, INTENT(inout) :: Npols
	567	INTEGER, DIMENSION(dx,dy), INTENT(inout) :: pols
	568	LOGICAL, INTENT(in) :: dbg
	569
	570	! Local
	571	INTEGER :: i,j,ip,iprm
	572	INTEGER, DIMENSION(Npols) :: origPol, NotPol, neigPol
	573	INTEGER :: ispol, NnotPol
	574	CHARACTER(len=4) :: ISa
	575
	576	!!!!!!! Variables
	577	! dx, dy: size of the space
	578	! Lmat: original bolean matrix from which the polygons come from
	579	! Npols: original number of polygons
	580	! pols: polygons space
	581
	582	fname = 'clean_polygons'
	583	IF (dbg) PRINT *," At '" // TRIM(fname) // "' ..."
	584
	585	origPol = -1
	586
	587	! Looking for polygons already in space
	588	NnotPol = 0
	589	DO ip=1, Npols
	590	ispol = COUNT(pols-ip == 0)
	591	IF (ispol > 0) THEN
	592	origPol(ip) = ip
	593	ELSE
	594	NnotPol = NnotPol + 1
	595	NotPol(NnotPol) = ip
	596	neigPol(NnotPol) = -1
	597	END IF
	598	END DO
	599
	600	IF (dbg) THEN
	601	PRINT *,' It should be:', Npols, ' polygons, but already there are:', Npols - NnotPol
	602	PRINT *,' Polygons to remove:', NotPol(1:NnotPol)
	603	END IF
	604
	605	! Looking for the hole border of a polygon. This is identify as such polygon point which along
	606	! y-axis has NpolygonA, Npolygon, .FALSE.
	607	DO i=1,dx
	608	DO j=2,dy-1
	609	IF ( (pols(i,j-1) /= pols(i,j) .AND. pols(i,j+1) == -1) .AND. (COUNT(NotPol-pols(i,j)==0)==0) &
	610	.AND. (pols(i,j) /= -1) .AND. (pols(i,j-1) /= -1)) THEN
	611	IF (dbg) PRINT *,' Polygon:', pols(i,j), ' to be removed at point (',i,',',j,'); j-1:', &
	612	pols(i,j-1), ' j:', pols(i,j), ' j+1:', pols(i,j+1)
	613	NnotPol = NnotPol + 1
	614	NotPol(NnotPol) = pols(i,j)
	615	neigPol(NnotPol) = pols(i,j-1)
	616	END IF
	617	END DO
	618	END DO
	619
	620	IF (dbg) THEN
	621	PRINT *,' It should be:', Npols, ' polygons, but already there are:', Npols - NnotPol
	622	PRINT *,' Polygons to remove after looking for fake border-of-hole polygons _______'
	623	DO i=1, NnotPol
	624	PRINT *, ' Polygon:', NotPol(i), ' to be replaced by:', neigPol(i)
	625	END DO
	626	END IF
	627
	628	! Removing polygons
	629	DO iprm=1, NnotPol
	630	IF (neigPol(iprm) == -1) THEN
	631	WHERE (pols == NotPol(iprm))
	632	pols = -1
	633	END WHERE
	634	IF (dbg) THEN
	635	PRINT *,' removing polygon:', NotPol(iprm)
	636	END IF
	637	ELSE
	638	WHERE (pols == NotPol(iprm))
	639	pols = neigPol(iprm)
	640	END WHERE
	641	IF (dbg) THEN
	642	PRINT *,' replacing polygon:', NotPol(iprm), ' by:', neigPol(iprm)
	643	END IF
	644	END IF
	645	END DO
	646
	647	! Re-numbering (descending values)
	648	DO i = 1, NnotPol
	649	iprm = MAXVAL(NotPol(1:NnotPol))
	650	WHERE(pols > iprm)
	651	pols = pols - 1
	652	END WHERE
	653	j = Index1DArrayI(NotPol, NnotPol, iprm)
	654	NotPol(j) = -9
	655	END DO
	656
	657	Npols = Npols - NnotPol
	658
	659	RETURN
	660
	661	END SUBROUTINE clean_polygons
	662
[1612]	663	SUBROUTINE path_properties(dx, dy, Lmat, Nptspth, pth, xxtrm, yxtrm, meanctr, axs, Nvrtx, vrtxs)
	664	! Subroutine to determine the properties of a path:
	665	! extremes: minimum and maximum of the path along x,y axes
	666	! meancenter: center from the mean of the coordinates of the paths locations
	667	! vertexs: path point, without neighbours along a given axis
	668
	669	IMPLICIT NONE
	670
	671	INTEGER, INTENT(in) :: dx, dy, Nptspth
	672	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat
	673	INTEGER, DIMENSION(Nptspth,2), INTENT(in) :: pth
	674	CHARACTER, INTENT(in) :: axs
	675	INTEGER, DIMENSION(2), INTENT(out) :: meanctr, xxtrm, yxtrm
	676	INTEGER, INTENT(out) :: Nvrtx
	677	INTEGER, DIMENSION(Nptspth,2), INTENT(out) :: vrtxs
	678
	679	! Local
	680	INTEGER :: i, ip, jp
	681	INTEGER :: neig1, neig2
	682
	683	!!!!!!! Variables
	684	! dx,dy: size of the space
	685	! Lmat: original matrix of logical values for the path
	686	! Nptspth: number of points of the path
	687	! pth: path coordinates (clockwise)
	688	! axs: axis of finding the vertex
	689	! [x/y]xtrm: minimum and maximum coordinates of the path
	690	! meanctr: center from the mean of the coordinates of the path
	691	! Nvrtx: Number of vertexs of the path
	692	! vrtxs: coordinates of the vertexs
	693
	694	fname = 'path_properties'
	695
	696	vrtxs = -1
	697	Nvrtx = 0
	698
	699	xxtrm = (/ MINVAL(pth(:,1)), MAXVAL(pth(:,1)) /)
	700	yxtrm = (/ MINVAL(pth(:,2)), MAXVAL(pth(:,2)) /)
	701	meanctr = (/ SUM(pth(:,1))/Nptspth, SUM(pth(:,2))/Nptspth /)
	702
	703	IF (axs == 'x' .OR. axs == 'X') THEN
	704	! Looking vertexs along x-axis
	705	DO i=1, Nptspth
	706	ip = pth(i,1)
	707	jp = pth(i,2)
	708	neig1 = 0
	709	neig2 = 0
	710	! W-point
	711	IF (ip == 1) THEN
	712	neig1 = -1
	713	ELSE
	714	IF (.NOT.Lmat(ip-1,jp)) neig1 = -1
	715	END IF
	716	! E-point
	717	IF (ip == dx) THEN
	718	neig2 = -1
	719	ELSE
	720	IF (.NOT.Lmat(ip+1,jp)) neig2 = -1
	721	END IF
	722
	723	IF (neig1 == -1 .AND. neig2 == -1) THEN
	724	Nvrtx = Nvrtx + 1
	725	vrtxs(Nvrtx,:) = (/ip,jp/)
	726	END IF
	727	END DO
	728	ELSE IF (axs == 'y' .OR. axs == 'Y') THEN
	729	! Looking vertexs along x-axis
	730	DO i=1, Nptspth
	731	ip = pth(i,1)
	732	jp = pth(i,2)
	733
	734	neig1 = 0
	735	neig2 = 0
	736	! S-point
	737	IF (jp == 1) THEN
	738	neig1 = -1
	739	ELSE
	740	IF (.NOT.Lmat(ip,jp-1)) neig1 = -1
	741	END IF
	742	! N-point
	743	IF (jp == dy) THEN
	744	neig2 = -1
	745	ELSE
	746	IF (.NOT.Lmat(ip,jp+1)) neig2 = -1
	747	END IF
	748
	749	IF (neig1 == -1 .AND. neig2 == -1) THEN
	750	Nvrtx = Nvrtx + 1
	751	vrtxs(Nvrtx,:) = (/ ip, jp /)
	752	END IF
	753	END DO
	754	ELSE
	755	msg = "Axis '" // axs // "' not available" // CHAR(10) // " Available ones: 'x', 'X', 'y, 'Y'"
	756	CALL ErrMsg(msg, fname, -1)
	757	END IF
	758
	759	RETURN
	760
	761	END SUBROUTINE path_properties
	762
	763	SUBROUTINE gridpoints_InsidePolygon(dx, dy, isbrdr, Npath, path, Nvrtx, xpathxtrm, ypathxtrm, &
	764	vrtxs, Npts, pts, inside)
	765	! Subroutine to determine if a series of grid points are inside a polygon following ray casting algorithm
	766	! FROM: https://en.wikipedia.org/wiki/Point_in_polygon
	767
	768	IMPLICIT NONE
	769
	770	INTEGER, INTENT(in) :: dx,dy,Npath,Nvrtx,Npts
	771	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isbrdr
	772	INTEGER, DIMENSION(Npath,2), INTENT(in) :: path
	773	INTEGER, DIMENSION(2), INTENT(in) :: xpathxtrm, ypathxtrm
	774	INTEGER, DIMENSION(Npath,2) :: vrtxs
	775	INTEGER, DIMENSION(Npts,2), INTENT(in) :: pts
	776	LOGICAL, DIMENSION(Npts), INTENT(out) :: inside
	777
	778	! Local
	779	INTEGER :: i,j,ip,ix,iy
[1614]	780	INTEGER :: Nintersecs, isvertex, ispath
[1612]	781	INTEGER :: ierr
	782	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: halo_brdr
[1614]	783	INTEGER :: Nbrbrdr
[1612]	784
	785	!!!!!!! Variables
	786	! dx,dy: space size
	787	! Npath: number of points of the path of the polygon
	788	! path: path of the polygon
	789	! isbrdr: boolean matrix of the space wqith .T. on polygon border
	790	! Nvrtx: number of vertexs of the path
[1614]	791	! [x/y]pathxtrm extremes of the path
[1612]	792	! vrtxs: vertexs of the path along y-axis
	793	! Npts: number of points
	794	! pts: points to look for
	795	! inside: vector wether point is inside or not (coincident to a border is inside)
	796
	797	fname = 'gridpoints_InsidePolygon'
	798
	799	! Creation of a 1-grid point larger matrix to deal with points reaching the limits
	800	IF (ALLOCATED(halo_brdr)) DEALLOCATE(halo_brdr)
	801	ALLOCATE(halo_brdr(dx+2,dy+2), STAT=ierr)
	802	msg = "Problems allocating matrix 'halo_brdr'"
	803	CALL ErrMsg(msg, fname, ierr)
	804	halo_brdr = .FALSE.
	805
	806	DO i=1,dx
	807	halo_brdr(i+1,2:dy+1) = isbrdr(i,:)
	808	END DO
	809
	810	inside = .FALSE.
	811
	812	DO ip=1,Npts
	813	Nintersecs = 0
	814	ix = pts(ip,1)
	815	iy = pts(ip,2)
	816	! Point might be outside path range...
	817	IF (ix >= xpathxtrm(1) .AND. ix <= xpathxtrm(2) .AND. iy >= ypathxtrm(1) .AND. &
	818	iy <= ypathxtrm(2)) THEN
	819
	820	! It is a border point?
[1614]	821	ispath = index_list_coordsI(Npath, path, (/ix,iy/))
	822	IF (isbrdr(ix,iy) .AND. (ispath /= -1)) THEN
[1612]	823	inside(ip) = .TRUE.
	824	CYCLE
	825	END IF
	826
	827	! Looking along y-axis
[1614]	828	! Accounting for consecutives borders
	829	Nbrbrdr = 0
	830	DO j=MAX(1,ypathxtrm(1)-1),iy-1
[1612]	831	! Only counting that borders that are not vertexs
[1614]	832	ispath = index_list_coordsI(Npath, path, (/ix,j/))
[1612]	833	isvertex = index_list_coordsI(Npath, vrtxs, (/ix,j/))
[1614]	834
	835	IF (halo_brdr(ix+1,j+1) .AND. (ispath /= -1) .AND. (isvertex == -1) ) Nintersecs = Nintersecs + 1
	836	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
	837	Nbrbrdr = Nbrbrdr + 1
	838	ELSE
	839	! Will remove that consecutive borders above 2
	840	IF (Nbrbrdr /= 0) THEN
	841	Nintersecs = Nintersecs - MAX(Nbrbrdr-1, 0)
	842	Nbrbrdr = 0
	843	END IF
	844	END IF
[1612]	845	END DO
	846	IF (MOD(Nintersecs,2) /= 0) inside(ip) = .TRUE.
	847	END IF
	848
	849	END DO
	850
	851	RETURN
	852
	853	END SUBROUTINE gridpoints_InsidePolygon
	854
	855	SUBROUTINE look_clockwise_borders(dx,dy,Nbrdrs,brdrs,gbrdr,isbrdr,ix,iy,dbg,xf,yf,iff)
	856	! Subroutine to look clock-wise for a next point within a collection of borders (limits of a region)
	857
	858	IMPLICIT NONE
	859
	860	INTEGER, INTENT(in) :: dx, dy, Nbrdrs, ix, iy
	861	INTEGER, DIMENSION(Nbrdrs,2), INTENT(in) :: brdrs
	862	LOGICAL, DIMENSION(Nbrdrs), INTENT(in) :: gbrdr
	863	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isbrdr
	864	LOGICAL, INTENT(in) :: dbg
	865	INTEGER, INTENT(out) :: xf, yf, iff
	866
	867	! Local
	868	INTEGER :: isch
	869	CHARACTER(len=2), DIMENSION(8) :: Lclock
	870	INTEGER, DIMENSION(8,2) :: spt
	871	INTEGER :: iif, jjf
	872
	873	!!!!!!! Variables
	874	! dx, dy: 2D shape ot the space
	875	! Nbrdrs: number of brdrs found in this 2D space
	876	! brdrs: list of coordinates of the borders
	877	! gbrdr: accounts for the use if the given border point
	878	! isbrdr: accounts for the matrix of the point is a border or not
	879	! ix,iy: coordinates of the point to start to find for
	880	! xf,yf: coordinates of the found point
	881	! iff: position of the border found within the list of borders
	882
	883	fname = 'look_clockwise_borders'
	884
	885	! Looking clock-wise assuming that one starts from the westernmost point
	886
	887	! Label of the search
[1614]	888	lclock = (/ 'W ', 'NW', 'N ', 'NE', 'E ', 'SE', 'S ', 'SW' /)
[1612]	889	! Transformation to apply
[1614]	890	!spt = (/ (/-1,0/), (/-1,1/), (/0,1/), (/1,1/), (/1,0/), (/1,-1/), (/0,-1/), (/-1,-1/) /)
	891	spt(:,1) = (/ -1, -1, 0, 1, 1, 1, 0, -1 /)
	892	spt(:,2) = (/ 0, 1, 1, 1, 0, -1, -1, -1 /)
[1612]	893
	894	xf = -1
	895	yf = -1
	896	DO isch=1, 8
	897	! clock-wise search
	898	IF (spt(isch,1) >= 0) THEN
	899	iif = MIN(dx,ix+spt(isch,1))
	900	ELSE
	901	iif = MAX(1,ix+spt(isch,1))
	902	END IF
	903	IF (spt(isch,2) >= 0) THEN
	904	jjf = MIN(dy,iy+spt(isch,2))
	905	ELSE
	906	jjf = MAX(1,iy+spt(isch,2))
	907	END IF
	908	iff = index_list_coordsI(Nbrdrs, brdrs,(/iif,jjf/))
	909	IF (iff > 0) THEN
	910	IF (dbg) PRINT *,' ' // lclock(isch) // '-point:', iif,jjf, ':', iff, 'is',isbrdr(iif,jjf), &
	911	'got',gbrdr(iff)
	912	IF (isbrdr(iif,jjf) .AND. .NOT.gbrdr(iff)) THEN
	913	xf = iif
	914	yf = jjf
	915	EXIT
	916	END IF
	917	END IF
	918	END DO
	919
	920	RETURN
	921
	922	END SUBROUTINE look_clockwise_borders
	923
[1614]	924	SUBROUTINE borders_matrixL(dx,dy,dxy,Lmat,brdrs,isbrdr,isbrdry)
[1612]	925	! Subroutine to provide the borders of a logical array (interested in .TRUE.)
	926
	927	IMPLICIT NONE
	928
	929	INTEGER, INTENT(in) :: dx,dy,dxy
	930	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat
	931	INTEGER, DIMENSION(dxy,2), INTENT(out) :: brdrs
[1614]	932	LOGICAL, DIMENSION(dx,dy), INTENT(out) :: isbrdr, isbrdry
[1612]	933
	934	! Local
	935	INTEGER :: i,j,ib
	936
	937	!!!!!!! Variables
	938	! dx,dy: size of the space
	939	! dxy: maximum number of border points
	940	! Lmat: Matrix to look for the borders
	941	! brdrs: list of coordinates of the borders
	942	! isbrdr: matrix with .T./.F. wether the given matrix point is a border or not
[1614]	943	! isbrdry: matrix with .T./.F. wether the given matrix point is a border or not only along y-axis
[1612]	944
	945	fname = 'borders_matrixL'
	946
	947	isbrdr = .FALSE.
	948	brdrs = -1
	949	ib = 1
	950
	951	! Starting with the borders. If a given point is TRUE it is a path-vertex
	952	! Along y-axis
	953	DO i=1, dx
	954	IF (Lmat(i,1) .AND. .NOT.isbrdr(i,1)) THEN
	955	brdrs(ib,1) = i
	956	brdrs(ib,2) = 1
	957	isbrdr(i,1) = .TRUE.
	958	ib=ib+1
	959	END IF
	960	IF (Lmat(i,dy) .AND. .NOT.isbrdr(i,dy)) THEN
	961	brdrs(ib,1) = i
	962	brdrs(ib,2) = dy
	963	isbrdr(i,dy) = .TRUE.
	964	ib=ib+1
	965	END IF
	966	END DO
	967	! Along x-axis
	968	DO j=1, dy
	969	IF (Lmat(1,j) .AND. .NOT.isbrdr(1,j)) THEN
	970	brdrs(ib,1) = 1
	971	brdrs(ib,2) = j
	972	isbrdr(1,j) = .TRUE.
	973	ib=ib+1
	974	END IF
	975	IF (Lmat(dx,j) .AND. .NOT.isbrdr(dx,j)) THEN
	976	brdrs(ib,1) = dx
	977	brdrs(ib,2) = j
	978	isbrdr(dx,j) = .TRUE.
	979	ib=ib+1
	980	END IF
	981	END DO
	982
[1614]	983	isbrdry = isbrdr
	984
[1612]	985	! Border as that when looking on x-axis points with Lmat(i) /= Lmat(i+1)
	986	DO i=1, dx-1
	987	DO j=1, dy-1
	988	IF ( Lmat(i,j) .NEQV. Lmat(i+1,j) ) THEN
	989	IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN
	990	brdrs(ib,1) = i
	991	brdrs(ib,2) = j
	992	isbrdr(i,j) = .TRUE.
	993	ib=ib+1
	994	ELSE IF (Lmat(i+1,j) .AND. .NOT.isbrdr(i+1,j)) THEN
	995	brdrs(ib,1) = i+1
	996	brdrs(ib,2) = j
	997	isbrdr(i+1,j) = .TRUE.
	998	ib=ib+1
	999	END IF
	1000	END IF
	1001	! y-axis
	1002	IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN
	1003	IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN
	1004	brdrs(ib,1) = i
	1005	brdrs(ib,2) = j
	1006	isbrdr(i,j) = .TRUE.
[1614]	1007	isbrdry(i,j) = .TRUE.
[1612]	1008	ib=ib+1
	1009	ELSE IF (Lmat(i,j+1) .AND. .NOT.isbrdr(i,j+1)) THEN
	1010	brdrs(ib,1) = i
	1011	brdrs(ib,2) = j+1
	1012	isbrdr(i,j+1) = .TRUE.
[1614]	1013	isbrdry(i,j+1) = .TRUE.
[1612]	1014	ib=ib+1
	1015	END IF
	1016	END IF
	1017	END DO
	1018	END DO
	1019
[1614]	1020	DO i=1, dx-1
	1021	DO j=1, dy-1
	1022	! y-axis
	1023	IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN
	1024	IF (Lmat(i,j)) THEN
	1025	isbrdry(i,j) = .TRUE.
	1026	ELSE IF (Lmat(i,j+1)) THEN
	1027	isbrdry(i,j+1) = .TRUE.
	1028	END IF
	1029	END IF
	1030	END DO
	1031	END DO
	1032	! only y-axis adding bands of 2 grid points
	1033	DO i=1, dx-1
	1034	DO j=2, dy-2
	1035	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
	1036	IF (Lmat(i,j)) THEN
	1037	isbrdry(i,j) = .TRUE.
	1038	isbrdry(i,j+1) = .TRUE.
	1039	END IF
	1040	END IF
	1041	END DO
	1042	END DO
	1043
[1612]	1044	RETURN
	1045
	1046	END SUBROUTINE borders_matrixL
	1047
	1048	SUBROUTINE paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)
	1049	! Subroutine to search the paths of a border field.
	1050
	1051	IMPLICIT NONE
	1052
	1053	INTEGER, INTENT(in) :: dx, dy, Nppt
	1054	LOGICAL, INTENT(in) :: dbg
	1055	LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isborder
	1056	INTEGER, DIMENSION(Nppt,2), INTENT(in) :: borders
	1057	INTEGER, DIMENSION(Nppt,Nppt,2), INTENT(out) :: paths
	1058	INTEGER, INTENT(out) :: Npath
	1059	INTEGER, DIMENSION(Nppt), INTENT(out) :: Nptpaths
	1060
	1061	! Local
	1062	INTEGER :: i,j,k,ib
	1063	INTEGER :: ierr
	1064	INTEGER :: Nbrdr
	1065	LOGICAL, DIMENSION(:), ALLOCATABLE :: gotbrdr, emptygotbrdr
	1066	INTEGER :: iipth, ipath, ip, Nptspath
	1067	INTEGER :: iib, jjb, iip, ijp, iif, jjf, iff
	1068	LOGICAL :: found, finishedstep
	1069
	1070	!!!!!!! Variables
	1071	! dx,dy: spatial dimensions of the space
	1072	! Nppt: possible number of paths and points that the paths can have
	1073	! isborder: boolean matrix which provide the borders of the polygon
	1074	! borders: coordinates of the borders of the polygon
	1075	! paths: coordinates of each found path
	1076	! Npath: number of paths found
	1077	! Nptpaths: number of points per path
	1078
	1079	fname = 'paths_border'
	1080
	1081	! Sarting matrix
	1082	paths = -1
	1083	Npath = 0
	1084	Nptspath = 0
	1085	Nptpaths = -1
	1086
	1087	ib=1
	1088	finishedstep = .FALSE.
	1089
	1090	! Number of border points
	1091	DO ib=1, Nppt
	1092	IF (borders(ib,1) == -1 ) EXIT
	1093	END DO
	1094	Nbrdr = ib-1
	1095
	1096	IF (dbg) THEN
	1097	PRINT *,' borders _______'
	1098	DO i=1,Nbrdr
	1099	PRINT *,' ',i,':',borders(i,:)
	1100	END DO
	1101	END IF
	1102
	1103	! Matrix which keeps track if a border point has been located
	1104	IF (ALLOCATED(gotbrdr)) DEALLOCATE(gotbrdr)
	1105	ALLOCATE(gotbrdr(Nbrdr), STAT=ierr)
	1106	msg = "Problems allocating matrix 'gotbrdr'"
	1107	CALL ErrMsg(msg, fname, ierr)
	1108	IF (ALLOCATED(emptygotbrdr)) DEALLOCATE(emptygotbrdr)
	1109	ALLOCATE(emptygotbrdr(Nbrdr), STAT=ierr)
	1110	msg = "Problems allocating matrix 'emptygotbrdr'"
	1111	CALL ErrMsg(msg, fname, ierr)
	1112
	1113	gotbrdr = .FALSE.
	1114	emptygotbrdr = .FALSE.
	1115
	1116	! Starting the fun...
	1117
	1118	! Looking along the lines and when a border is found, starting from there in a clock-wise way
	1119	iipth = 1
	1120	ipath = 1
	1121	DO ib=1,Nbrdr
	1122	iib = borders(iipth,1)
	1123	jjb = borders(iipth,2)
	1124	! Starting new path
	1125	newpath: IF (.NOT.gotbrdr(iipth)) THEN
	1126	ip = 1
	1127	Nptspath = 1
	1128	paths(ipath,ip,:) = borders(iipth,:)
	1129	gotbrdr(iipth) = .TRUE.
	1130	! Looking for following clock-wise search
	1131	! Not looking for W, because search starts from the W
	1132	iip = iib
	1133	ijp = jjb
	1134	DO k=1,Nbrdr
	1135	IF (dbg) PRINT *,ipath,'iip jip:', iip, ijp
	1136	found = .FALSE.
	1137	CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE.,iif,jjf,iff)
	1138	IF (iif /= -1) THEN
	1139	ip=ip+1
	1140	paths(ipath,ip,:) = (/ iif,jjf /)
	1141	found = .TRUE.
	1142	gotbrdr(iff) = .TRUE.
	1143	iip = iif
	1144	ijp = jjf
	1145	Nptspath = Nptspath + 1
	1146	END IF
	1147
	1148	IF (dbg) THEN
	1149	PRINT *,iib,jjb,' end of this round path:', ipath, '_____', gotbrdr
	1150	DO i=1, Nptspath
	1151	PRINT *,' ',i,':',paths(ipath,i,:)
	1152	END DO
	1153	END IF
	1154	! If it is not found a next point, might be because it is a non-polygon related value
	1155	IF (.NOT.found) THEN
	1156	IF (dbg) PRINT *,'NOT FOUND !!!', gotbrdr
	1157	! Are still there available borders?
	1158	IF (ALL(gotbrdr) .EQV. .TRUE.) THEN
	1159	finishedstep = .TRUE.
	1160	Npath = ipath
	1161	Nptpaths(ipath) = Nptspath
	1162	EXIT
	1163	ELSE
	1164	Nptpaths(ipath) = Nptspath
[1614]	1165	! Let's have a look if the previous points in the path have already some 'non-located' neighbourgs
	1166	DO i=Nptspath,1,-1
	1167	iip = paths(ipath,i,1)
	1168	ijp = paths(ipath,i,2)
	1169	CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE., iif, &
	1170	jjf,iff)
	1171	IF (iif /= -1 .AND. iff /= -1) THEN
	1172	IF (dbg) PRINT *,' re-take path from point:', iif,',',jjf,' n-path:', iff
	1173	found = .TRUE.
	1174	iipth = index_list_coordsI(Nppt, borders, (/iip,ijp/))
	1175	EXIT
[1612]	1176	END IF
[1614]	1177	END DO
	1178	IF (.NOT.found) THEN
[1612]	1179	! Looking for the next available border point for the new path
	1180	DO i=1,Nbrdr
	1181	IF (.NOT.gotbrdr(i)) THEN
	1182	iipth = i
	1183	EXIT
	1184	END IF
	1185	END DO
	1186	IF (dbg) PRINT *,' Looking for next path starting at:', iipth, ' point:', &
	1187	borders(iipth,:)
	1188	ipath=ipath+1
	1189	EXIT
	1190	END IF
	1191	END IF
	1192	ELSE
	1193	IF (dbg) PRINT *,' looking for next point...'
	1194	END IF
	1195	IF (finishedstep) EXIT
	1196	END DO
	1197	END IF newpath
	1198	END DO
	1199	Npath = ipath
	1200	Nptpaths(ipath) = Nptspath
	1201
	1202	DEALLOCATE (gotbrdr)
	1203	DEALLOCATE (emptygotbrdr)
	1204
	1205	RETURN
	1206
	1207	END SUBROUTINE paths_border
	1208
[1609]	1209	SUBROUTINE rand_sample(Nvals, Nsample, sample)
	1210	! Subroutine to randomly sample a range of indices
	1211
	1212	IMPLICIT NONE
	1213
	1214	INTEGER, INTENT(in) :: Nvals, Nsample
	1215	INTEGER, DIMENSION(Nsample), INTENT(out) :: sample
	1216
	1217	! Local
	1218	INTEGER :: i, ind, jmax
	1219	REAL, DIMENSION(Nsample) :: randv
	1220	CHARACTER(len=50) :: fname
	1221	LOGICAL :: found
	1222	LOGICAL, DIMENSION(Nvals) :: issampled
	1223	CHARACTER(len=256) :: msg
	1224	CHARACTER(len=10) :: IS1, IS2
	1225
	1226	!!!!!!! Variables
	1227	! Nvals: number of values
	1228	! Nsamples: number of samples
	1229	! sample: samnple
	1230	fname = 'rand_sample'
	1231
	1232	IF (Nsample > Nvals) THEN
	1233	WRITE(IS1,'(I10)')Nvals
	1234	WRITE(IS2,'(I10)')Nsample
	1235	msg = 'Sampling of ' // TRIM(IS1) // ' is too big for ' // TRIM(IS1) // 'values'
	1236	CALL ErrMsg(msg, fname, -1)
	1237	END IF
	1238
	1239	! Generation of random numbers always the same series during the whole program!
	1240	CALL RANDOM_NUMBER(randv)
	1241
	1242	! Making sure that we do not repeat any value
	1243	issampled = .FALSE.
	1244
	1245	DO i=1, Nsample
	1246	! Generation of the index from the random numbers
	1247	ind = MAX(INT(randv(i)*Nvals), 1)
	1248
	1249	IF (.NOT.issampled(ind)) THEN
	1250	sample(i) = ind
	1251	issampled(ind) = .TRUE.
	1252	ELSE
	1253	! Looking around the given index
	1254	!PRINT *,' Index :', ind, ' already sampled!', issampled(ind)
	1255	found = .FALSE.
	1256	DO jmax=1, Nvals
	1257	ind = MIN(ind+jmax, Nvals)
	1258	IF (.NOT.issampled(ind)) THEN
	1259	sample(i) = ind
	1260	issampled(ind) = .TRUE.
	1261	found = .TRUE.
	1262	EXIT
	1263	END IF
	1264	ind = MAX(1, ind-jmax)
	1265	IF (.NOT.issampled(ind)) THEN
	1266	sample(i) = ind
	1267	issampled(ind) = .TRUE.
	1268	found = .TRUE.
	1269	EXIT
	1270	END IF
	1271	END DO
	1272	IF (.NOT.found) THEN
	1273	msg = 'sampling could not be finished due to absence of available value!!'
	1274	CALL ErrMsg(msg, fname, -1)
	1275	END IF
	1276	END IF
	1277
	1278	END DO
	1279
	1280	RETURN
	1281
	1282	END SUBROUTINE rand_sample
	1283
	1284	SUBROUTINE PrintQuantilesR_K(Nvals, vals, Nquants, qtvs, bspc)
	1285	! Subroutine to print the quantiles of values REAL(r_k)
	1286
	1287	IMPLICIT NONE
	1288
	1289	INTEGER, INTENT(in) :: Nvals, Nquants
	1290	REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals
	1291	REAL(r_k), DIMENSION(Nquants), INTENT(in) :: qtvs
	1292	CHARACTER(len=1000), OPTIONAL :: bspc
	1293
	1294	! Local
	1295	INTEGER :: iq
	1296	LOGICAL, DIMENSION(Nvals) :: search1, search2, search
	1297	CHARACTER(len=6) :: RS1
	1298	CHARACTER(len=50) :: fname
	1299	CHARACTER(len=1000) :: bspcS
	1300
	1301	!!!!!!! Variables
	1302	! vals: series of values
	1303	! qtvs: values of the quantiles
	1304	! bspc: base quantity of spaces
	1305
	1306	fname = 'PrintQuantilesR_K'
	1307
	1308	IF (PRESENT(bspc)) THEN
	1309	bspcS = bspc
	1310	ELSE
	1311	bspcS = ' '
	1312	END IF
	1313
	1314	DO iq=1, Nquants-1
	1315
	1316	WHERE (vals >= qtvs(iq))
	1317	search1 = .TRUE.
	1318	ELSEWHERE
	1319	search1 = .FALSE.
	1320	END WHERE
	1321
	1322	WHERE (vals < qtvs(iq+1))
	1323	search2 = .TRUE.
	1324	ELSEWHERE
	1325	search2 = .FALSE.
	1326	END WHERE
	1327
	1328	WHERE (search1 .AND. search2)
	1329	search = .TRUE.
	1330	ELSEWHERE
	1331	search = .FALSE.
	1332	END WHERE
	1333
	1334	WRITE(RS1, '(F6.2)')(iq)*100./(Nquants-1)
	1335	PRINT *, TRIM(bspcS) // '[',iq,']', TRIM(RS1) // ' %:', qtvs(iq), 'N:', COUNT(search)
	1336
	1337	END DO
	1338
	1339	RETURN
	1340
	1341	END SUBROUTINE PrintQuantilesR_K
	1342
	1343	INTEGER FUNCTION FindMinimumR_K(x, dsize, Startv, Endv)
	1344	! Function returns the location of the minimum in the section between Start and End.
	1345
	1346	IMPLICIT NONE
	1347
	1348	INTEGER, INTENT(in) :: dsize
	1349	REAL(r_k), DIMENSION(dsize), INTENT(in) :: x
	1350	INTEGER, INTENT(in) :: Startv, Endv
	1351
	1352	! Local
	1353	REAL(r_k) :: Minimum
	1354	INTEGER :: Location
	1355	INTEGER :: i
	1356
	1357	Minimum = x(Startv) ! assume the first is the min
	1358	Location = Startv ! record its position
	1359	DO i = Startv+1, Endv ! start with next elements
	1360	IF (x(i) < Minimum) THEN ! if x(i) less than the min?
	1361	Minimum = x(i) ! Yes, a new minimum found
	1362	Location = i ! record its position
	1363	END IF
	1364	END DO
	1365
	1366	FindMinimumR_K = Location ! return the position
	1367
	1368	END FUNCTION FindMinimumR_K
	1369
	1370	SUBROUTINE SwapR_K(a, b)
	1371	! Subroutine swaps the values of its two formal arguments.
	1372
	1373	IMPLICIT NONE
	1374
	1375	REAL(r_k), INTENT(INOUT) :: a, b
	1376	! Local
	1377	REAL(r_k) :: Temp
	1378
	1379	Temp = a
	1380	a = b
	1381	b = Temp
	1382
	1383	END SUBROUTINE SwapR_K
	1384
	1385	SUBROUTINE SortR_K(x, Nx)
	1386	! Subroutine receives an array x() r_K and sorts it into ascending order.
	1387
	1388	IMPLICIT NONE
	1389
	1390	INTEGER, INTENT(IN) :: Nx
	1391	REAL(r_k), DIMENSION(Nx), INTENT(INOUT) :: x
	1392
	1393	! Local
	1394	INTEGER :: i
	1395	INTEGER :: Location
	1396
	1397	DO i = 1, Nx-1 ! except for the last
	1398	Location = FindMinimumR_K(x, Nx-i+1, i, Nx) ! find min from this to last
	1399	CALL SwapR_K(x(i), x(Location)) ! swap this and the minimum
	1400	END DO
	1401
	1402	END SUBROUTINE SortR_K
	1403
	1404	SUBROUTINE quantilesR_K(Nvals, vals, Nquants, quants)
	1405	! Subroutine to provide the quantiles of a given set of values of type real 'r_k'
	1406
	1407	IMPLICIT NONE
	1408
	1409	INTEGER, INTENT(in) :: Nvals, Nquants
	1410	REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals
	1411	REAL(r_k), DIMENSION(Nquants), INTENT(out) :: quants
	1412
	1413	! Local
	1414	INTEGER :: i
	1415	REAL(r_k) :: minv, maxv
	1416	REAL(r_k), DIMENSION(Nvals) :: sortedvals
	1417
	1418	!!!!!!! Variables
	1419	! Nvals: number of values
	1420	! Rk: kind of real of the values
	1421	! vals: values
	1422	! Nquants: number of quants
	1423	! quants: values at which the quantile start
	1424
	1425	minv = MINVAL(vals)
	1426	maxv = MAXVAL(vals)
	1427
	1428	sortedvals = vals
	1429	! Using from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
	1430	CALL SortR_K(sortedvals, Nvals)
	1431
	1432	quants(1) = minv
	1433	DO i=2, Nquants
	1434	quants(i) = sortedvals(INT((i-1)*Nvals/Nquants))
	1435	END DO
	1436
	1437	END SUBROUTINE quantilesR_K
	1438
[1618]	1439
	1440	SUBROUTINE StatsR_K(Nvals, vals, minv, maxv, mean, mean2, stdev)
	1441	! Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a
	1442	! series of r_k numbers
	1443
	1444	IMPLICIT NONE
	1445
	1446	INTEGER, INTENT(in) :: Nvals
	1447	REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals
	1448	REAL(r_k), INTENT(out) :: minv, maxv, mean, mean2, stdev
	1449
	1450	!!!!!!! Variables
	1451	! Nvals: number of values
	1452	! vals: values
	1453	! minv: minimum value of values
	1454	! maxv: maximum value of values
	1455	! mean: mean value of values
	1456	! mean2: quadratic mean value of values
	1457	! stdev: standard deviation of values
	1458
	1459	minv = MINVAL(vals)
	1460	maxv = MAXVAL(vals)
	1461
	1462	mean=SUM(vals)
	1463	mean2=SUM(vals*vals)
	1464
	1465	mean=mean/Nvals
	1466	mean2=mean2/Nvals
	1467
	1468	stdev=SQRT(mean2 - mean*mean)
	1469
	1470	RETURN
	1471
	1472	END SUBROUTINE StatsR_k
	1473
[1609]	1474	END MODULE module_scientific

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

Download in other formats:

Original Format