67 integer,
pointer :: ipt(:) => null(), ipl(:) => null()
68 integer :: nt=imiss,nl=imiss
72 double precision :: x,y
78 MODULE PROCEDURE triangles_delete
81INTERFACE triangles_compute
82 MODULE PROCEDURE triangles_compute_r, triangles_compute_d, triangles_compute_c
86public triangles, triangles_new,
delete, triangles_compute, xy
92function triangles_new(ndp)
result(this)
93type(triangles) :: this
94integer,
intent(in) :: ndp
101if (
c_e(ndp) .and. ndp >= 3)
then
102 allocate(this%ipt(6*ndp-15), this%ipl(6*ndp))
108end function triangles_new
112subroutine triangles_delete(this)
113type(triangles) :: this
115if (
associated(this%ipt))
deallocate(this%ipt)
116if (
associated(this%ipl))
deallocate(this%ipl)
121end subroutine triangles_delete
124integer function triangles_compute_r (XD,YD,tri)
125real,
intent(in) :: XD(:)
126real,
intent(in) :: YD(:)
127type (triangles),
intent(inout) :: tri
128type (xy) :: co(size(xd))
133 triangles_compute_r = contng_simc(co,tri%NT,tri%IPT,tri%NL,tri%IPL)
135end function triangles_compute_r
137integer function triangles_compute_d (XD,YD,tri)
138double precision,
intent(in) :: XD(:)
139double precision,
intent(in) :: YD(:)
140type (triangles),
intent(inout) :: tri
141type (xy) :: co(size(xd))
146 triangles_compute_d = contng_simc(co,tri%NT,tri%IPT,tri%NL,tri%IPL)
148end function triangles_compute_d
150integer function triangles_compute_c (co,tri)
151type (xy),
intent(in) :: co(:)
152type (triangles),
intent(inout) :: tri
155 triangles_compute_c = contng_simc(co,tri%NT,tri%IPT,tri%NL,tri%IPL)
157end function triangles_compute_c
173integer function contng_simc (co,NT,IPT,NL,IPL)
175type (xy),
intent(in) :: co(:)
176integer,
intent(out) :: NT
177integer,
intent(out) :: NL
182integer,
intent(out) :: IPT(:)
188integer,
intent(out) :: IPL(:)
197integer :: IWL(18*size(co)), IWP(size(co))
199double precision :: WK(size(co))
200integer :: ITF(2), NREP=100
202double precision :: AR,ARMN,ARMX,DSQ12,DSQI,DSQMN,DSQMX,DX,DX21,DXMN,DXMX,DY,DY21,DYMN,DYMX
203double precision :: X1,Y1
204integer :: ilf,IP,IP1,IP1P1,IP2,IP3,IPL1,IPL2,IPLJ1,IPLJ2,IPMN1,IPMN2,IPT1,IPT2,IPT3
205INTEGER :: IPTI,IPTI1,IPTI2,IREP,IT,IT1T3,IT2T3,ITS,ITT3
206INTEGER :: ILFT2,ITT3R,JLT3,JP,JP1,JP2,JP2T3,JP3T3,JPC,JPMN,JPMX,JWL,JWL1
207INTEGER :: JWL1MN,NDP,NDPM1,NL0,NLF,NLFC,NLFT2,NLN,NLNT3,NLT3,NSH,NSHT3,NT0,NTF
208INTEGER :: NTT3,NTT3P3
211integer ::i,mloc(size(co)-1),mlocall(1),mlocv(1)
228call l4f_log(l4f_debug,
"start triangulation")
231 mlocv=minloc(vdsqf(co(i),co(i+1:)))+i
235mlocall=minloc((/(vdsqf(co(i),co(mloc(i))),i=1,
size(mloc))/))
237dsqmn = vdsqf(co(mlocall(1)),co(mloc(mlocall(1))))
239ipmn2 = mloc(mlocall(1))
241call l4f_log(l4f_debug,
"end triangulation closest pair")
248 call l4f_log(l4f_error,
"CONTNG-IDENTICAL INPUT DATA POINTS FOUND")
285dmp%x = (co(ipmn1)%x+co(ipmn2)%x)/2.0
286dmp%y = (co(ipmn1)%y+co(ipmn2)%y)/2.0
294 IF (ip1.EQ.ipmn1 .OR. ip1.EQ.ipmn2) cycle
297 wk(jp1) = vdsqf(dmp,co(ip1))
303 IF (wk(jp2) .GE. dsqmn) cycle
313call l4f_log(l4f_debug,
"end triangulation sort")
328 IF (abs((co(ip)%y-y1)*dx21-(co(ip)%x-x1)*dy21) .GT. ar)
then
334 call l4f_log(l4f_debug,
"CONTNG - ALL COLLINEAR DATA POINTS")
347call l4f_log(l4f_debug,
"end triangulation collinear")
358IF (side(co(ip1),co(ip2),co(ip3)) < 10.0)
then
380call l4f_log(l4f_debug,
"end triangulation first triangle")
396 dsqmn = dxmn**2+dymn**2
410 IF (ar.GE.(-armn) .AND. dsqi.GE.dsqmn)
GO TO 230
417230 ar = dy*dxmx-dx*dymx
418 IF (ar .LT. (-armx)) cycle
420 IF (ar.LE.armx .AND. dsqi.GE.dsqmx) cycle
427 IF (jpmx .LT. jpmn) jpmx = jpmx+nl0
437 ipl(jp3t3-2) = ipl(jp2t3-2)
438 ipl(jp3t3-1) = ipl(jp2t3-1)
439 ipl(jp3t3) = ipl(jp2t3)
443 ipl(jp2t3-2) = ipl(jp3t3-2)
444 ipl(jp2t3-1) = ipl(jp3t3-1)
445 ipl(jp2t3) = ipl(jp3t3)
456 l310 :
DO jp2=jpmx,nl0
472 IF (jp2 == jpmx)
then
480 ipl(nlnt3-1) = ipl(1)
489 IF (ipti.NE.ipl1 .AND. ipti.NE.ipl2)
GO TO 300
491 IF (ipti.NE.ipl1 .AND. ipti.NE.ipl2)
GO TO 300
496300
IF (conxch_simc(co%X,co%Y,ip1,ipti,ipl1,ipl2) .EQ. 0) cycle l310
504 IF (jp2 .EQ. jpmx) ipl(jp2t3) = it
505 IF (jp2.EQ.nl0 .AND. ipl(3).EQ.it) ipl(3) = nt0
518 IF (nlf .EQ. 0) cycle l400
538 IF (ipl1.NE.ipt1 .AND. ipl1.NE.ipt2 .AND. ipl1.NE.ipt3) cycle
539 IF (ipl2.NE.ipt1 .AND. ipl2.NE.ipt2 .AND. ipl2.NE.ipt3) cycle
542 IF (ntf .EQ. 2)
GO TO 330
544 IF (ntf .LT. 2) cycle l370
551 IF (ipti1.NE.ipl1 .AND. ipti1.NE.ipl2)
GO TO 340
553 IF (ipti1.NE.ipl1 .AND. ipti1.NE.ipl2)
GO TO 340
557 IF (ipti2.NE.ipl1 .AND. ipti2.NE.ipl2)
GO TO 350
559 IF (ipti2.NE.ipl1 .AND. ipti2.NE.ipl2)
GO TO 350
564350
IF (conxch_simc(co%X,co%Y,ipti1,ipti2,ipl1,ipl2) .EQ. 0) cycle l370
589 IF ((iplj1.EQ.ipl1 .AND. iplj2.EQ.ipti2) .OR. &
590 (iplj2.EQ.ipl1 .AND. iplj1.EQ.ipti2)) ipl(jlt3) = itf(1)
591 IF ((iplj1.EQ.ipl2 .AND. iplj2.EQ.ipti1) .OR. &
592 (iplj2.EQ.ipl2 .AND. iplj1.EQ.ipti1)) ipl(jlt3) = itf(2)
597 IF (nlf .EQ. nlfc) cycle l400
604 DO jwl1=jwl1mn,nlft2,2
606 iwl(jwl-1) = iwl(jwl1-1)
613call l4f_log(l4f_debug,
"end triangulation appending")
623 IF (side(co(ip1),co(ip2),co(ip3)) .GE. 10.0) cycle
627call l4f_log(l4f_debug,
"end triangulation rearranging")
632call l4f_log(l4f_debug,
"end triangulation")
645elemental double precision function vdsqf(co1,co2)
646type(xy),
intent(in) :: co1,co2
648vdsqf = (co2%x-co1%x)**2+(co2%y-co1%y)**2
649if (vdsqf == 0.d0) vdsqf = huge(vdsqf)
659double precision function side(co1,co2,co3)
660type(xy),
intent(in):: co1,co2,co3
662side = (co3%y-co1%y)*(co2%x-co1%x)-(co3%x-co1%x)*(co2%y-co1%y)
665end function contng_simc
673INTEGER FUNCTION conxch_simc (X,Y,I1,I2,I3,I4)
675double precision,
intent(in) :: x(:), y(:)
678integer,
intent(in) :: i1,i2,i3,i4
680double precision :: x0(4), y0(4)
681double precision :: c2sq,c1sq,a3sq,b2sq,b3sq,a1sq,a4sq,b1sq,b4sq,a2sq,c4sq,c3sq
683double precision :: s1sq,s2sq,s3sq,s4sq,u1,u2,u3,u4
685equivalence(c2sq,c1sq),(a3sq,b2sq),(b3sq,a1sq),(a4sq,b1sq), &
686 (b4sq,a2sq),(c4sq,c3sq)
697u3 = (y0(2)-y0(3))*(x0(1)-x0(3))-(x0(2)-x0(3))*(y0(1)-y0(3))
698u4 = (y0(1)-y0(4))*(x0(2)-x0(4))-(x0(1)-x0(4))*(y0(2)-y0(4))
700 u1 = (y0(3)-y0(1))*(x0(4)-x0(1))-(x0(3)-x0(1))*(y0(4)-y0(1))
701 u2 = (y0(4)-y0(2))*(x0(3)-x0(2))-(x0(4)-x0(2))*(y0(3)-y0(2))
702 a1sq = (x0(1)-x0(3))**2+(y0(1)-y0(3))**2
703 b1sq = (x0(4)-x0(1))**2+(y0(4)-y0(1))**2
704 c1sq = (x0(3)-x0(4))**2+(y0(3)-y0(4))**2
705 a2sq = (x0(2)-x0(4))**2+(y0(2)-y0(4))**2
706 b2sq = (x0(3)-x0(2))**2+(y0(3)-y0(2))**2
707 c3sq = (x0(2)-x0(1))**2+(y0(2)-y0(1))**2
708 s1sq = u1*u1/(c1sq*max(a1sq,b1sq))
709 s2sq = u2*u2/(c2sq*max(a2sq,b2sq))
710 s3sq = u3*u3/(c3sq*max(a3sq,b3sq))
711 s4sq = u4*u4/(c4sq*max(a4sq,b4sq))
712 IF (min(s1sq,s2sq) < min(s3sq,s4sq)) idx = 1
717END FUNCTION conxch_simc
Function to check whether a value is missing or not.
Distructor for triangles.
Utilities for CHARACTER variables.
classe per la gestione del logging
Definitions of constants and functions for working with missing values.
Space utilities, derived from NCAR software.