Context Navigation

module_scientific.f90 @ 1626

Last change on this file since 1626 was 1622, checked in by lfita, 8 years ago
Fixing various issues
File size: 49.3 KB

Line
1	MODULE module_scientific
2	! Module of the scientific function/subroutines
3
4	!!!!!!! Functions & subroutines
5	! all_polygons_properties: Subroutine to determine the properties of all polygons in a 2D field:
6	! borders_matrixL: Subroutine to provide the borders of a logical array (interested in .TRUE.)
7	! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
8	! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
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
15	! polygon_properties: Subroutine to determine the properties of a polygon (as .TRUE. matrix)
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!
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.
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
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
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)
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
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
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
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
51
52	IMPLICIT NONE
53
54	LOGICAL, INTENT(in) :: dbg
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
59	CHARACTER(len=1000), INTENT(in) :: projN
60	INTEGER, DIMENSION(Npoly), INTENT(out) :: Npolyptss
61	INTEGER, DIMENSION(Npoly,2), INTENT(out) :: xxtrms, yxtrms, meanctrs
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
65	INTEGER, DIMENSION(Npoly,2), INTENT(out) :: nvcoords, xvcoords
66	REAL(r_k), DIMENSION(Npoly,2), INTENT(out) :: meanwctrs, meanvnctrs, meanvxctrs
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
81	! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
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
85	! meanwctrs: lon, lat coordinates of the center from the spatial-weighted mean of the polygons
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
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
95
96	fname = 'all_polygons_properties'
97
98	! Initializing matrices
99	Npolyptss = -1
100	xxtrms = fillval64
101	yxtrms = fillval64
102	meanctrs = fillval64
103	meanwctrs = fillval64
104	areas = fillval64
105	nvals = fillvalI
106	xvals = fillval64
107	mvals = fillval64
108	m2vals = fillval64
109	stdvals = fillval64
110	quants = fillval64
111	nvcoords = fillvalI
112	xvcoords = fillvalI
113	meanvnctrs = fillval64
114	meanvxctrs = fillval64
115
116	DO ip=1, Npoly
117	boolpoly = polys == ip
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,:))
122	END DO
123
124	RETURN
125
126	END SUBROUTINE all_polygons_properties
127
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)
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
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
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
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
145
146	IMPLICIT NONE
147
148	LOGICAL, INTENT(in) :: dbg
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
153	CHARACTER(len=1000), INTENT(in) :: projN
154	INTEGER, INTENT(out) :: Npolypts
155	INTEGER, DIMENSION(2), INTENT(out) :: xxtrm, yxtrm, meanctr
156	INTEGER, DIMENSION(2), INTENT(out) :: nvcoord, xvcoord
157	REAL(r_k), DIMENSION(2), INTENT(out) :: meanwctr, meanvnctr, meanvxctr
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
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
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
182	! 'ctsarea': Constant Area
183	! 'lon/lat': for regular longitude-latitude
184	! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
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
188	! meanwctr: lon, lat coordinates of the center from the spatial-weighted mean of the polygon
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
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
198
199	fname = 'polygon_properties'
200
201	IF (dbg) PRINT *," '" // TRIM(fname) // "' ..."
202
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
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
253	ip = 1
254	DO i=1,dx
255	DO j=1,dy
256	IF (poly(i,j)) THEN
257	polypts(ip,1) = i
258	polypts(ip,2) = j
259	polyvals(ip) = values(i,j)
260	SELECT CASE (TRIM(projS))
261	CASE ('ctsarea')
262	! Constant Area
263	xcorr(i,j) = 1.
264	ycorr(i,j) = 1.
265	CASE ('lon/lat')
266	! Area as fixed yres and sinus-lat varying for xres
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.
274	END SELECT
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
288	ip = ip + 1
289	END IF
290	END DO
291	END DO
292
293	IF (dbg) THEN
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
299	END IF
300
301	sumRESx = xres*SUM(xcorr)
302	sumRESy = yres*SUM(ycorr)
303
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 /)
307	meanwctr = (/ meanwctr(1)/sumRESx, meanwctr(2)/sumRESy /)
308
309	IF (dbg) THEN
310	PRINT *,' mean grid center: ', meanctr, ' weighted mean center: ', meanwctr
311	END IF
312
313	! Statistics of the values within the polygon
314	CALL StatsR_K(Npolypts, polyvals, nval, xval, mval, m2val, stdval)
315
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
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)
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
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
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
400	Npoly = 0
401
402	DO it=1,dt
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
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
435	LOGICAL, DIMENSION(dx,dy) :: isborder, isbordery
436	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: paths
437	INTEGER :: Npath
438	INTEGER, DIMENSION(:), ALLOCATABLE :: Nptpaths
439	INTEGER, DIMENSION(2) :: xtrx, xtry, meanpth
440	INTEGER :: Nvertx, Npts
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
454	! The mathematical maximum woiuld be dx*dy/4, but let's be optimistic... (sorry Jero)
455	Nppt = dx*dy/10
456
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
472	! Filling with the points of all the space with .TRUE.
473	Npts = COUNT(boolmat)
474
475	IF (ALLOCATED(points)) DEALLOCATE(points)
476	ALLOCATE(points(Npts,2), STAT=ierr)
477	msg = "Problems allocating matrix 'points'"
478	CALL ErrMsg(msg, fname, ierr)
479
480	! We only want to localize that points 'inside'
481	ip = 1
482	DO i=1, dx
483	DO j=1, dy
484	IF (boolmat(i,j)) THEN
485	points(ip,1) = i
486	points(ip,2) = j
487	ip = ip + 1
488	END IF
489	END DO
490	END DO
491
492	CALL borders_matrixL(dx, dy, Nppt, boolmat, borders, isborder, isbordery)
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)
504	ALLOCATE(isin(Npts), STAT=ierr)
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
526	CALL gridpoints_InsidePolygon(dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), Nvertx,&
527	xtrx, xtry, vertxs, Npts, points, isin)
528
529	! Filling polygons
530	DO ipp=1, Npts
531	IF (isin(ipp)) polys(points(ipp,1),points(ipp,2)) = ip
532	END DO
533
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
543	END DO
544
545	! Cleaning polygons matrix of not-used and paths around holes
546	CALL clean_polygons(dx, dy, boolmat, polys, Npoly, dbg)
547
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
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
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
780	INTEGER :: Nintersecs, isvertex, ispath
781	INTEGER :: ierr
782	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: halo_brdr
783	INTEGER :: Nbrbrdr
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
791	! [x/y]pathxtrm extremes of the path
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?
821	ispath = index_list_coordsI(Npath, path, (/ix,iy/))
822	IF (isbrdr(ix,iy) .AND. (ispath /= -1)) THEN
823	inside(ip) = .TRUE.
824	CYCLE
825	END IF
826
827	! Looking along y-axis
828	! Accounting for consecutives borders
829	Nbrbrdr = 0
830	DO j=MAX(1,ypathxtrm(1)-1),iy-1
831	! Only counting that borders that are not vertexs
832	ispath = index_list_coordsI(Npath, path, (/ix,j/))
833	isvertex = index_list_coordsI(Npath, vrtxs, (/ix,j/))
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
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
888	lclock = (/ 'W ', 'NW', 'N ', 'NE', 'E ', 'SE', 'S ', 'SW' /)
889	! Transformation to apply
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 /)
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
924	SUBROUTINE borders_matrixL(dx,dy,dxy,Lmat,brdrs,isbrdr,isbrdry)
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
932	LOGICAL, DIMENSION(dx,dy), INTENT(out) :: isbrdr, isbrdry
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
943	! isbrdry: matrix with .T./.F. wether the given matrix point is a border or not only along y-axis
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
983	isbrdry = isbrdr
984
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.
1007	isbrdry(i,j) = .TRUE.
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.
1013	isbrdry(i,j+1) = .TRUE.
1014	ib=ib+1
1015	END IF
1016	END IF
1017	END DO
1018	END DO
1019
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
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
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
1176	END IF
1177	END DO
1178	IF (.NOT.found) THEN
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
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
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
1474	END MODULE module_scientific

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

Download in other formats:

Original Format