Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

polygones.f90 @ 5143

Last change on this file since 5143 was 3965, checked in by jan.polcher, 8 years ago

Merge with trunk at revision3959.
This includes all the developments made for CMIP6 and passage to XIOS2.
All conflicts are resolved and the code compiles.

But it still does not link because of an "undefined reference to `_intel_fast_memmove'"

File size: 22.9 KB

Line
1
2	MODULE polygones
3	!
4	! This module is there to handle polygones. It is an essential tool in interpolation methods.
5	!
6	USE defprec
7
8	IMPLICIT NONE
9
10	PRIVATE
11	PUBLIC :: polygones_pointinside, polygones_extend, polygones_intersection, polygones_cleanup, &
12	& polygones_area, polygones_crossing, polygones_convexhull
13	!
14	REAL(r_std),PARAMETER :: zero=0.0
15	!
16	CONTAINS
17	!
18	!=======================================================================================================
19	!
20	SUBROUTINE polygones_pointinside(nvert_in, poly, point_x, point_y, inside)
21	!
22	! Routine to determine of point (point_x, point_y) is inside of polygone poly
23	! Mathode based on docume http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
24	!
25	! Arguments
26	INTEGER(i_std), INTENT(in) :: nvert_in
27	REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly
28	REAL(r_std), INTENT(in) :: point_x, point_y
29	LOGICAL, INTENT(out) :: inside
30	!
31	! Local
32	INTEGER(i_std) :: i, j, nvert
33	REAL(r_std) :: xonline
34	!
35	!
36	inside = .FALSE.
37	!
38	nvert=SIZE(poly,DIM=1)
39	IF ( nvert < nvert_in) THEN
40	WRITE(,) "Input polygone is smaller than the size proposed in the interface."
41	STOP "polygones_pointinside"
42	ENDIF
43	IF (SIZE(poly,DIM=2) .NE. 2) THEN
44	WRITE(,) "This cannot be a polygone :", SIZE(poly)
45	STOP "polygones_pointinside"
46	ENDIF
47	!
48	j = nvert_in
49	DO i=1,nvert_in
50	IF ( (poly(i,2) > point_y) .NEQV. (poly(j,2) > point_y) ) THEN
51	xonline = (poly(j,1)-poly(i,1))*(point_y-poly(i,2))/(poly(j,2)-poly(i,2))+poly(i,1)
52	IF ( point_x < xonline ) THEN
53	inside = (.NOT. inside)
54	ENDIF
55	ENDIF
56	j = i
57	ENDDO
58	!
59	END SUBROUTINE polygones_pointinside
60	!
61	!=======================================================================================================
62	!
63	SUBROUTINE polygones_lineintersect(la, lb, intersection, point_x, point_y)
64	!
65	! Simple alogorith tho determine the intersections of 2 lines.
66	!
67	! ARGUMENTS
68	REAL(r_std), DIMENSION(2,2), INTENT(in) :: la, lb
69	LOGICAL, INTENT(out) :: intersection
70	REAL(r_std), INTENT(out) :: point_x, point_y
71	!
72	! Local
73	REAL(r_std) :: den, ua, ub
74	REAL(r_std) :: x1, x2, x3, x4, y1, y2, y3, y4
75	!
76	intersection = .FALSE.
77	!
78	x1 = la(1,1)
79	y1 = la(1,2)
80	x2 = la(2,1)
81	y2 = la(2,2)
82
83	x3 = lb(1,1)
84	y3 = lb(1,2)
85	x4 = lb(2,1)
86	y4 = lb(2,2)
87
88	den = (y4-y3)(x2-x1)-(x4-x3)(y2-y1)
89
90	IF ( ABS(den) > EPSILON(den) ) THEN
91	ua = ((x4-x3)(y1-y3)-(y4-y3)(x1-x3))/den
92	ub = ((x2-x1)(y1-y3)-(y2-y1)(x1-x3))/den
93	IF ( ua >= 0 .AND. ua <= 1 .AND. ub >= 0 .AND. ub <= 1) THEN
94	intersection = .TRUE.
95	point_x = x1 + ua*(x2-x1)
96	point_y = y1 + ua*(y2-y1)
97	ENDIF
98	ENDIF
99
100	END SUBROUTINE polygones_lineintersect
101	!
102	!=======================================================================================================
103	!
104	SUBROUTINE polygones_extend(nvert_in, poly_in, nbdots, nvert_out, poly_out)
105	!
106	! This function simply add nbdots onto each side of the polygone.
107	!
108	! Arguments
109	INTEGER(i_std), INTENT(in) :: nvert_in
110	REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_in
111	INTEGER(i_std), INTENT(in) :: nbdots
112	INTEGER(i_std), INTENT(out) :: nvert_out
113	REAL(r_std), DIMENSION(:,:), INTENT(out) :: poly_out
114	!
115	! Local
116	!
117	INTEGER(i_std) :: nvert, nvert_tmp, i, j, id, ipos
118	REAL(r_std) :: xs, xe, ys, ye
119	!
120	nvert=SIZE(poly_in,DIM=1)
121	IF ( nvert < nvert_in ) THEN
122	WRITE(,) "Input polygone is smaller than the size proposed in the interface."
123	STOP "polygones_extend"
124	ENDIF
125	IF (SIZE(poly_in,DIM=2) .NE. 2) THEN
126	WRITE(,) "Poly_in cannot be a polygone :", SIZE(poly_in)
127	STOP "polygones_extend"
128	ENDIF
129	nvert_tmp=SIZE(poly_out,DIM=1)
130	IF (SIZE(poly_out,DIM=2) .NE. 2) THEN
131	WRITE(,) "Poly_out cannot be a polygone :", SIZE(poly_out)
132	STOP "polygones_extend"
133	ENDIF
134	!
135	IF ( nvert_tmp < nvert_in*nbdots) THEN
136	WRITE(,) "Poly_out is too small for extending it :", SIZE(poly_out), " we would need :", nvert_in*nbdots
137	STOP "polygones_extend"
138	ENDIF
139	nvert_out = nvert_in*nbdots
140	!
141	j = nvert_in
142	ipos = 0
143	DO i=1,nvert_in
144	xs = poly_in(j,1)
145	xe = poly_in(i,1)
146	ys = poly_in(j,2)
147	ye = poly_in(i,2)
148	IF (xs == xe ) THEN
149	DO id=1,nbdots
150	ipos = ipos + 1
151	poly_out(ipos,1) = xs
152	poly_out(ipos,2) = ys + (id-1)*(ye-ys)/nbdots
153	ENDDO
154	ELSE IF (ys == ye ) THEN
155	DO id=1,nbdots
156	ipos = ipos + 1
157	poly_out(ipos,1) = xs + (id-1)*(xe-xs)/nbdots
158	poly_out(ipos,2) = ys
159	ENDDO
160	ELSE
161	DO id=1,nbdots
162	ipos = ipos + 1
163	poly_out(ipos,2) = ys + (id-1)*(ye-ys)/nbdots
164	poly_out(ipos,1) = (xs-xe)*(poly_out(ipos,2)-ye)/(ys-ye)+xe
165	ENDDO
166	ENDIF
167	j = i
168	ENDDO
169	!
170	END SUBROUTINE polygones_extend
171	!
172	!=======================================================================================================
173	!
174	SUBROUTINE polygones_intersection(nvert_a, poly_a, nvert_b, poly_b, nbpts, poly_out)
175	!
176	! Find the polygon resulting from the intersection of poly_a and poly_b
177	!
178	! Arguments
179	INTEGER(i_std), INTENT(in) :: nvert_a
180	REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_a
181	INTEGER(i_std), INTENT(in) :: nvert_b
182	REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_b
183	INTEGER(i_std), INTENT(out) :: nbpts
184	REAL(r_std), DIMENSION(:,:), INTENT(out) :: poly_out
185	!
186	! Local
187	INTEGER(i_std) :: nvert_tmpa, nvert_tmpb, nvert_out, i, j
188	LOGICAL :: inside
189	INTEGER(i_std) :: in_a, in_b
190	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: ptin_a, ptin_b
191	!
192	nvert_tmpa=SIZE(poly_a,DIM=1)
193	IF ( nvert_tmpa < nvert_a) THEN
194	WRITE(,) "Input polygone (poly_a) is smaller than the size proposed in the interface."
195	STOP "polygones_intersection"
196	ENDIF
197	IF (SIZE(poly_a,DIM=2) .NE. 2) THEN
198	WRITE(,) "Poly_a cannot be a polygone :", SIZE(poly_a)
199	STOP "polygones_intersection"
200	ENDIF
201	ALLOCATE(ptin_a(nvert_a))
202	nvert_tmpb=SIZE(poly_b,DIM=1)
203	IF ( nvert_tmpb < nvert_a) THEN
204	WRITE(,) "Input polygone (poly_b) is smaller than the size proposed in the interface."
205	STOP "polygones_intersection"
206	ENDIF
207	IF (SIZE(poly_b,DIM=2) .NE. 2) THEN
208	WRITE(,) "Poly_b cannot be a polygone :", SIZE(poly_b)
209	STOP "polygones_intersection"
210	ENDIF
211	ALLOCATE(ptin_b(nvert_b))
212	!
213	!
214	in_a = 0
215	DO i=1,nvert_a
216	CALL polygones_pointinside(nvert_b, poly_b, poly_a(i,1), poly_a(i,2), inside)
217	IF ( inside ) THEN
218	in_a = in_a + 1
219	ptin_a(in_a) = i
220	ENDIF
221	ENDDO
222	!
223	in_b = 0
224	DO i=1,nvert_b
225	CALL polygones_pointinside(nvert_a, poly_a, poly_b(i,1), poly_b(i,2), inside)
226	IF ( inside ) THEN
227	in_b = in_b + 1
228	ptin_b(in_b) = i
229	ENDIF
230	ENDDO
231	!
232	nbpts = in_a + in_b
233	!
234	nvert_out=SIZE(poly_out,DIM=1)
235	IF (SIZE(poly_out,DIM=2) .NE. 2) THEN
236	WRITE(,) "Poly_out cannot be a polygone :", SIZE(poly_out)
237	STOP "polygones_intersection"
238	ENDIF
239	IF ( nbpts <= nvert_out ) THEN
240	i = 0
241	DO j=1,in_a
242	i = i+1
243	poly_out(i,:) = poly_a(ptin_a(j),:)
244	ENDDO
245	DO j=1,in_b
246	i = i+1
247	poly_out(i,:) = poly_b(ptin_b(j),:)
248	ENDDO
249	ELSE
250	WRITE(,) "The intersection polygone has", nbpts, "points and that is more than available in poly_out", SIZE(poly_out)
251	STOP "polygones_intersection"
252	ENDIF
253	!
254	DEALLOCATE(ptin_a, ptin_b)
255	!
256	END SUBROUTINE polygones_intersection
257	!
258	!=======================================================================================================
259	!
260	SUBROUTINE polygones_cleanup(nvert_in, poly_in, nvert_out, poly_out)
261	!
262	! Clean_up the polygone by deleting points which are redundant and ordering based
263	! on the proximity of the points. Beware this does not give a convex hull to the surface
264	! inside the polygone.
265	!
266	! Arguments
267	INTEGER(i_std), INTENT(in) :: nvert_in
268	REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_in
269	INTEGER(i_std), INTENT(out) :: nvert_out
270	REAL(r_std), DIMENSION(:,:), INTENT(out) :: poly_out
271	!
272	! Local
273	!
274	INTEGER(i_std) :: nvert, nvert_tmpin, nvert_tmp
275	INTEGER(i_std) :: i, j
276	INTEGER(i_std), DIMENSION(1) :: ismin
277	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: dist
278	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: poly_tmp
279	!
280	nvert_tmpin=SIZE(poly_in,DIM=1)
281	IF ( nvert_tmpin < nvert_in ) THEN
282	WRITE(,) "Input polygone is smaller than the size proposed in the interface."
283	STOP "polygones_cleanup"
284	ENDIF
285	IF (SIZE(poly_in,DIM=2) .NE. 2) THEN
286	WRITE(,) "Poly_in cannot be a polygone :", SIZE(poly_in)
287	STOP "polygones_cleanup"
288	ENDIF
289	nvert_tmp=SIZE(poly_out,DIM=1)
290	IF (SIZE(poly_out,DIM=2) .NE. 2) THEN
291	WRITE(,) "Poly_out cannot be a polygone :", SIZE(poly_out)
292	STOP "polygones_cleanup"
293	ENDIF
294	!
295	ALLOCATE(poly_tmp(nvert_in,2))
296	ALLOCATE(dist(nvert_in))
297	!
298	nvert = nvert_in-1
299	poly_tmp(1:nvert,:) = poly_in(2:nvert_in,:)
300	poly_out(1,:) = poly_in(1,:)
301	nvert_out = 1
302	!
303	DO WHILE ( nvert > 0 )
304	!
305	DO j=1,nvert
306	dist(j) = SQRT((poly_out(nvert_out,1)-poly_tmp(j,1))2+(poly_out(nvert_out,2)-poly_tmp(j,2))2)
307	ENDDO
308	!
309	ismin = MINLOC(dist(1:nvert))
310	!
311	IF ( dist(ismin(1)) > EPSILON(dist(1)) ) THEN
312	!
313	! The smallest distance between 2 vertices is larger than zero :
314	! Add the vertex to poly_out and delete it from poly_tmp
315	!
316	nvert_out = nvert_out + 1
317	IF ( nvert_out > nvert_tmp) THEN
318	WRITE(,) "Output polygone too small"
319	STOP "polygones_cleanup"
320	ENDIF
321	poly_out(nvert_out,:) = poly_tmp(ismin(1),:)
322	poly_tmp(ismin(1):nvert-1,:) = poly_tmp(ismin(1)+1:nvert,:)
323	ELSE
324	!
325	! Else the vertex already exists in poly_out and thus we just need
326	! to delete it from poly_tmp
327	!
328	poly_tmp(ismin(1):nvert-1,:) = poly_tmp(ismin(1)+1:nvert,:)
329	ENDIF
330	!
331	! One less vertices exists in poly_tmp
332	!
333	nvert = nvert-1
334	!
335	ENDDO
336
337	DEALLOCATE(poly_tmp, dist)
338
339	END SUBROUTINE polygones_cleanup
340	!
341	!=======================================================================================================
342	!
343	SUBROUTINE polygones_area(nvert_in, poly_in, dx, dy, area)
344	!
345	! Find the polygon resulting from the intersection of poly_a and poly_b
346	!
347	! Arguments
348	INTEGER(i_std), INTENT(in) :: nvert_in
349	REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_in
350	REAL(r_std), INTENT(in) :: dx, dy
351	REAL(r_std), INTENT(out) :: area
352	!
353	! Local
354	!
355	INTEGER(i_std) :: nvert, i, j
356	!
357	nvert = SIZE(poly_in,DIM=1)
358	IF ( nvert < nvert_in) THEN
359	WRITE(,) "Input polygone is smaller than the size proposed in the interface."
360	STOP "polygones_area"
361	ENDIF
362	IF (SIZE(poly_in,DIM=2) .NE. 2) THEN
363	WRITE(,) "Poly_in cannot be a polygone :", SIZE(poly_in)
364	STOP "polygones_area"
365	ENDIF
366	!
367	area = zero
368	j = nvert_in
369	DO i=1,nvert_in
370	!! area = area + (dy(poly_in(j,2)+poly_in(i,2))/2.0)(dx*(poly_in(j,1)-poly_in(i,1)))
371	area = area + dydx/2.0(poly_in(j,2)+poly_in(i,2))*(poly_in(j,1)-poly_in(i,1))
372	j = i
373	ENDDO
374	!
375	area = ABS(area)
376	!
377	END SUBROUTINE polygones_area
378	!
379	!=======================================================================================================
380	!
381	SUBROUTINE polygones_crossing(nvert_a, poly_a, nvert_b, poly_b, nbpts, crossings)
382	!
383	! Find the polygon resulting from the intersection of poly_a and poly_b
384	!
385	! Arguments
386	INTEGER(i_std), INTENT(in) :: nvert_a
387	REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_a
388	INTEGER(i_std), INTENT(in) :: nvert_b
389	REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_b
390	INTEGER(i_std), INTENT(out) :: nbpts
391	REAL(r_std), DIMENSION(:,:), INTENT(out) :: crossings
392	!
393	! Local
394	!
395	INTEGER(i_std) :: nvert, ia, ja, ib, jb
396	REAL(r_std), DIMENSION(2,2) :: la, lb
397	REAL(r_std) :: x, y
398	LOGICAL :: intersect
399	!
400	nvert=SIZE(poly_a,DIM=1)
401	IF ( nvert < nvert_a) THEN
402	WRITE(,) "Input polygone (poly_a) is smaller than the size proposed in the interface."
403	STOP "polygones_crossing"
404	ENDIF
405	IF (SIZE(poly_a,DIM=2) .NE. 2) THEN
406	WRITE(,) "Poly_a cannot be a polygone :", SIZE(poly_a)
407	STOP "polygones_crossing"
408	ENDIF
409	nvert=SIZE(poly_b,DIM=1)
410	IF ( nvert < nvert_b) THEN
411	WRITE(,) "Input polygone (poly_b) is smaller than the size proposed in the interface."
412	STOP "polygones_crossing"
413	ENDIF
414	IF (SIZE(poly_b,DIM=2) .NE. 2) THEN
415	WRITE(,) "Poly_b cannot be a polygone :", SIZE(poly_b)
416	STOP "polygones_crossing"
417	ENDIF
418	nvert=SIZE(crossings,DIM=1)
419	IF (SIZE(crossings,DIM=2) .NE. 2) THEN
420	WRITE(,) "Poly_a cannot be a polygone :", SIZE(poly_a)
421	STOP "polygones_crossing"
422	ENDIF
423	!
424	nbpts = 0
425	!
426	ja = nvert_a
427	DO ia=1,nvert_a
428	!
429	! Put one side of the polygone a into la
430	!
431	la(1,1) = poly_a(ja,1)
432	la(1,2) = poly_a(ja,2)
433	la(2,1) = poly_a(ia,1)
434	la(2,2) = poly_a(ia,2)
435	!
436	jb = nvert_b
437	DO ib=1,nvert_b
438	!
439	! Put one side of the polygone b into lb
440	!
441	lb(1,1) = poly_b(jb,1)
442	lb(1,2) = poly_b(jb,2)
443	lb(2,1) = poly_b(ib,1)
444	lb(2,2) = poly_b(ib,2)
445	!
446	CALL polygones_lineintersect(la, lb, intersect, x, y)
447	!
448	IF ( intersect ) THEN
449	nbpts = nbpts + 1
450	IF ( nbpts <= nvert ) THEN
451	crossings(nbpts, 1) = x
452	crossings(nbpts, 2) = y
453	ELSE
454	WRITE(,) "Polygone to write the intersection points is too small", nvert, nbpts
455	STOP "polygones_crossing"
456	ENDIF
457	ENDIF
458	!
459	jb = ib
460	ENDDO
461	ja = ia
462	ENDDO
463	!
464
465	END SUBROUTINE polygones_crossing
466	!
467	!=======================================================================================================
468	!
469	SUBROUTINE polygones_convexhull(nvert_in, poly_in, nvert_out, poly_out)
470	!
471	! Routine orders points in poly_in so that they form a convex shape. This is based on the Graham scan
472	! algorith as implemented by Alan Miller : http://jblevins.org/mirror/amiller/
473	!
474	! The output polygone can be smaller than the input one as we are computing the convex envelope.
475	!
476	! Arguments
477	INTEGER(i_std), INTENT(in) :: nvert_in
478	REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly_in
479	INTEGER(i_std), INTENT(out) :: nvert_out
480	REAL(r_std), DIMENSION(:,:), INTENT(out) :: poly_out
481	!
482	! Local
483	!
484	INTEGER(i_std) :: nvert, nvert_tmp
485	INTEGER(i_std), ALLOCATABLE, DIMENSION(:):: iwk, next, vertex
486	REAL(r_std) :: xmin, temp, dy, dx2, dy1, dy2
487	REAL(r_std) :: x1, x2, y1, y2
488	REAL(r_std) :: xmax, ymin, ymax
489	REAL(r_std) :: dx, dx1
490	REAL(r_std) :: dmax, dmax1, dmax2
491	REAL(r_std) :: dist, dmin
492	INTEGER(i_std) :: i, i1, i2, j, jp1, jp2, i2save, i3, i2next
493	LOGICAL :: points_todo
494	INTEGER(i_std), PARAMETER :: nextinc=20
495	!
496	nvert_tmp=SIZE(poly_in,DIM=1)
497	IF ( nvert_tmp < nvert_in) THEN
498	WRITE(,) "Input polygone is smaller than the size proposed in the interface."
499	STOP "polygones_convexhull"
500	ENDIF
501	IF ( nvert_tmp <= 2 ) THEN
502	WRITE(,) "The input polygone is too small to compute a convex hull."
503	STOP "polygones_convexhull"
504	ENDIF
505	IF (SIZE(poly_in,DIM=2) .NE. 2) THEN
506	WRITE(,) "Poly_a cannot be a polygone :", SIZE(poly_in)
507	STOP "polygones_convexhull"
508	ENDIF
509	IF (SIZE(poly_out,DIM=2) .NE. 2) THEN
510	WRITE(,) "Poly_b cannot be a polygone :", SIZE(poly_out)
511	STOP "polygones_convexhull"
512	ENDIF
513	!
514	ALLOCATE(vertex(nvert_in))
515	ALLOCATE(iwk(nvert_in))
516	ALLOCATE(next(nvert_in*nextinc))
517	!
518	!
519	! Choose the points with smallest & largest x- values as the
520	! first two vertices of the polygon.
521	!
522	IF (poly_in(1,1) > poly_in(nvert_in,1)) THEN
523	vertex(1) = nvert_in
524	vertex(2) = 1
525	xmin = poly_in(nvert_in,1)
526	xmax = poly_in(1,1)
527	ELSE
528	vertex(1) = 1
529	vertex(2) = nvert_in
530	xmin = poly_in(1,1)
531	xmax = poly_in(nvert_in,1)
532	END IF
533
534	DO i = 2, nvert_in-1
535	temp = poly_in(i,1)
536	IF (temp < xmin) THEN
537	vertex(1) = i
538	xmin = temp
539	ELSE IF (temp > xmax) THEN
540	vertex(2) = i
541	xmax = temp
542	END IF
543	END DO
544	!
545	! Special case, xmax = xmin.
546	!
547	IF (xmax == xmin) THEN
548	IF (poly_in(1,2) > poly_in(nvert_in,2)) THEN
549	vertex(1) = nvert_in
550	vertex(2) = 1
551	ymin = poly_in(nvert_in,2)
552	ymax = poly_in(1,2)
553	ELSE
554	vertex(1) = 1
555	vertex(2) = nvert_in
556	ymin = poly_in(1,2)
557	ymax = poly_in(nvert_in,2)
558	END IF
559
560	DO i = 2, nvert_in-1
561	temp = poly_in(i,2)
562	IF (temp < ymin) THEN
563	vertex(1) = i
564	ymin = temp
565	ELSE IF (temp > ymax) THEN
566	vertex(2) = i
567	ymax = temp
568	END IF
569	END DO
570
571	nvert = 2
572	IF (ymax == ymin) nvert = 1
573	RETURN
574	END IF
575	!
576	! Set up two initial lists of points; those points above & those below the
577	! line joining the first two vertices. next(i) will hold the pointer to the
578	! point furthest from the line joining vertex(i) to vertex(i+1) on the left
579	! hand side.
580	!
581	i1 = vertex(1)
582	i2 = vertex(2)
583	iwk(i1) = -1
584	iwk(i2) = -1
585	dx = xmax - xmin
586	y1 = poly_in(i1,2)
587	dy = poly_in(i2,2) - y1
588	dmax = zero
589	dmin = zero
590	next(1) = -1
591	next(2) = -1
592	!
593	DO i = 1, nvert_in
594	IF (i == vertex(1) .OR. i == vertex(2)) CYCLE
595	dist = (poly_in(i,2) - y1)dx - (poly_in(i,1) - xmin)dy
596	IF (dist > zero) THEN
597	iwk(i1) = i
598	i1 = i
599	IF (dist > dmax) THEN
600	next(1) = i
601	dmax = dist
602	END IF
603	ELSE IF (dist < zero) THEN
604	iwk(i2) = i
605	i2 = i
606	IF (dist < dmin) THEN
607	next(2) = i
608	dmin = dist
609	END IF
610	END IF
611	END DO
612	!
613	! Ends of lists are indicated by pointers to -ve positions.
614	!
615	iwk(i1) = -1
616	iwk(i2) = -1
617	nvert = 2
618	!
619	j = 1
620	!
621	! Start of main process.
622	!
623	! Introduce new vertex between vertices j & j+1, if one has been found.
624	! Otherwise increase j. Exit if no more vertices.
625	!
626	!
627	points_todo = .TRUE.
628	DO WHILE ( points_todo )
629	!
630	DO WHILE ( points_todo .AND. next(j) < 0 )
631	IF (j == nvert) points_todo = .FALSE.
632	j = j + 1
633	ENDDO
634	!
635	IF ( points_todo ) THEN
636	!
637	jp1 = j + 1
638	IF ( jp1 >= nvert_in*nextinc) THEN
639	STOP "polygones_convexhull : please increase nextinc"
640	ENDIF
641	!
642	DO i = nvert, jp1, -1
643	vertex(i+1) = vertex(i)
644	next(i+1) = next(i)
645	END DO
646	jp2 = jp1 + 1
647	nvert = nvert + 1
648	IF (jp2 > nvert) jp2 = 1
649	i1 = vertex(j)
650	i2 = next(j)
651	i3 = vertex(jp2)
652	vertex(jp1) = i2
653	!
654	! Process the list of points associated with vertex j. New list at vertex j
655	! consists of those points to the left of the line joining it to the new
656	! vertex (j+1). Similarly for the list at the new vertex.
657	! Points on or to the right of these lines are dropped.
658	!
659	x1 = poly_in(i1,1)
660	x2 = poly_in(i2,1)
661	y1 = poly_in(i1,2)
662	y2 = poly_in(i2,2)
663	!
664	dx1 = x2 - x1
665	dx2 = poly_in(i3,1) - x2
666	dy1 = y2 - y1
667	dy2 = poly_in(i3,2) - y2
668	dmax1 = zero
669	dmax2 = zero
670	next(j) = -1
671	next(jp1) = -1
672	i2save = i2
673	i2next = iwk(i2)
674	i = iwk(i1)
675	iwk(i1) = -1
676	iwk(i2) = -1
677	!
678	DO WHILE ( i > 0 )
679	IF (i /= i2save) THEN
680	dist = (poly_in(i,2) - y1)dx1 - (poly_in(i,1) - x1)dy1
681	IF (dist > zero) THEN
682	iwk(i1) = i
683	i1 = i
684	IF (dist > DMAX1) THEN
685	next(j) = i
686	dmax1 = dist
687	END IF
688	ELSE
689	dist = (poly_in(i,2) - y2)dx2 - (poly_in(i,1) - x2)dy2
690	IF (dist > zero) THEN
691	iwk(i2) = i
692	i2 = i
693	IF (dist > dmax2) THEN
694	next(jp1) = i
695	dmax2 = dist
696	END IF
697	END IF
698	END IF
699	i = iwk(i)
700	ELSE
701	i = i2next
702	END IF
703	!
704	! Get next point from old list at vertex j.
705	!
706	ENDDO
707	!
708	! End lists with -ve values.
709	!
710	iwk(i1) = -1
711	iwk(i2) = -1
712	!
713	ENDIF
714	ENDDO
715	!
716	! Copy the polygone in he right order to the output variable
717	!
718	nvert_out = nvert
719	nvert_tmp=SIZE(poly_out,DIM=1)
720	IF ( nvert_tmp < nvert_out) THEN
721	WRITE(,) "Ouptput polygone is smaller than the size proposed in the interface."
722	STOP "polygones_convexhull"
723	ENDIF
724	!
725	DO i=1,nvert
726	poly_out(i,:) = poly_in(vertex(i),:)
727	ENDDO
728	!
729	DEALLOCATE(vertex, iwk, next)
730	!
731	END SUBROUTINE polygones_convexhull
732
733	END MODULE polygones

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_global/polygones.f90 @ 5143

Download in other formats: