67 integer,
pointer :: ipt(:) => null(), ipl(:) => null()
68 integer :: nt=imiss,nl=imiss
72 double precision :: x,y
78 MODULE PROCEDURE triangles_delete
81 INTERFACE triangles_compute
82 MODULE PROCEDURE triangles_compute_r, triangles_compute_d, triangles_compute_c
86 public triangles, triangles_new,
delete, triangles_compute, xy
92 function triangles_new(ndp)
result(this)
93 type(triangles) :: this
94 integer,
intent(in) :: ndp
101 if (
c_e(ndp) .and. ndp >= 3)
then 102 allocate(this%ipt(6*ndp-15), this%ipl(6*ndp))
108 end function triangles_new
112 subroutine triangles_delete(this)
113 type(triangles) :: this
115 if (
associated(this%ipt))
deallocate(this%ipt)
116 if (
associated(this%ipl))
deallocate(this%ipl)
121 end subroutine triangles_delete
124 integer function triangles_compute_r (XD,YD,tri)
125 real,
intent(in) :: XD(:)
126 real,
intent(in) :: YD(:)
127 type(triangles),
intent(inout) :: tri
128 type(xy) :: co(
size(xd))
130 if (tri%nt /= 0)
then 133 triangles_compute_r = contng_simc(co,tri%NT,tri%IPT,tri%NL,tri%IPL)
135 end function triangles_compute_r
137 integer function triangles_compute_d (XD,YD,tri)
138 double precision,
intent(in) :: XD(:)
139 double precision,
intent(in) :: YD(:)
140 type(triangles),
intent(inout) :: tri
141 type(xy) :: co(
size(xd))
143 if (tri%nt /= 0)
then 146 triangles_compute_d = contng_simc(co,tri%NT,tri%IPT,tri%NL,tri%IPL)
148 end function triangles_compute_d
150 integer function triangles_compute_c (co,tri)
151 type(xy),
intent(in) :: co(:)
152 type(triangles),
intent(inout) :: tri
154 if (tri%nt /= 0)
then 155 triangles_compute_c = contng_simc(co,tri%NT,tri%IPT,tri%NL,tri%IPL)
157 end function triangles_compute_c
173 integer function contng_simc (co,NT,IPT,NL,IPL)
175 type(xy),
intent(in) :: co(:)
176 integer,
intent(out) :: NT
177 integer,
intent(out) :: NL
182 integer,
intent(out) :: IPT(:)
188 integer,
intent(out) :: IPL(:)
197 integer :: IWL(18*size(co)), IWP(size(co))
199 double precision :: WK(size(co))
200 integer :: ITF(2), NREP=100
202 double precision :: AR,ARMN,ARMX,DSQ12,DSQI,DSQMN,DSQMX,DX,DX21,DXMN,DXMX,DY,DY21,DYMN,DYMX
203 double precision :: X1,Y1
204 integer :: ilf,IP,IP1,IP1P1,IP2,IP3,IPL1,IPL2,IPLJ1,IPLJ2,IPMN1,IPMN2,IPT1,IPT2,IPT3
205 INTEGER :: IPTI,IPTI1,IPTI2,IREP,IT,IT1T3,IT2T3,ITS,ITT3
206 INTEGER :: ILFT2,ITT3R,JLT3,JP,JP1,JP2,JP2T3,JP3T3,JPC,JPMN,JPMX,JWL,JWL1
207 INTEGER :: JWL1MN,NDP,NDPM1,NL0,NLF,NLFC,NLFT2,NLN,NLNT3,NLT3,NSH,NSHT3,NT0,NTF
208 INTEGER :: NTT3,NTT3P3
211 integer ::i,mloc(size(co)-1),mlocall(1),mlocv(1)
228 call l4f_log(l4f_debug,
"start triangulation")
231 mlocv=minloc(vdsqf(co(i),co(i+1:)))+i
235 mlocall=minloc((/(vdsqf(co(i),co(mloc(i))),i=1,
size(mloc))/))
237 dsqmn = vdsqf(co(mlocall(1)),co(mloc(mlocall(1))))
239 ipmn2 = mloc(mlocall(1))
241 call l4f_log(l4f_debug,
"end triangulation closest pair")
244 IF (dsqmn == 0.)
then 248 call l4f_log(l4f_error,
"CONTNG-IDENTICAL INPUT DATA POINTS FOUND")
285 dmp%x = (co(ipmn1)%x+co(ipmn2)%x)/2.0
286 dmp%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
313 call l4f_log(l4f_debug,
"end triangulation sort")
322 dx21 = co(ipmn2)%x-x1
323 dy21 = co(ipmn2)%y-y1
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")
347 call l4f_log(l4f_debug,
"end triangulation collinear")
358 IF (side(co(ip1),co(ip2),co(ip3)) < 10.0)
then 380 call 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
417 230 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
496 300
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
564 350
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)
613 call l4f_log(l4f_debug,
"end triangulation appending")
623 IF (side(co(ip1),co(ip2),co(ip3)) .GE. 10.0) cycle
627 call l4f_log(l4f_debug,
"end triangulation rearranging")
632 call l4f_log(l4f_debug,
"end triangulation")
645 elemental double precision function vdsqf(co1,co2)
646 type(xy),
intent(in) :: co1,co2
648 vdsqf = (co2%x-co1%x)**2+(co2%y-co1%y)**2
649 if (vdsqf == 0.d0) vdsqf = huge(vdsqf)
659 double precision function side(co1,co2,co3)
660 type(xy),
intent(in):: co1,co2,co3
662 side = (co3%y-co1%y)*(co2%x-co1%x)-(co3%x-co1%x)*(co2%y-co1%y)
665 end function contng_simc
673 INTEGER FUNCTION conxch_simc (X,Y,I1,I2,I3,I4)
675 double precision,
intent(in) :: x(:), y(:)
678 integer,
intent(in) :: i1,i2,i3,i4
680 double precision :: x0(4), y0(4)
681 double precision :: c2sq,c1sq,a3sq,b2sq,b3sq,a1sq,a4sq,b1sq,b4sq,a2sq,c4sq,c3sq
683 double precision :: s1sq,s2sq,s3sq,s4sq,u1,u2,u3,u4
685 equivalence(c2sq,c1sq),(a3sq,b2sq),(b3sq,a1sq),(a4sq,b1sq), &
686 (b4sq,a2sq),(c4sq,c3sq)
697 u3 = (y0(2)-y0(3))*(x0(1)-x0(3))-(x0(2)-x0(3))*(y0(1)-y0(3))
698 u4 = (y0(1)-y0(4))*(x0(2)-x0(4))-(x0(1)-x0(4))*(y0(2)-y0(4))
699 IF (u3*u4 > 0.0)
then 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
717 END FUNCTION conxch_simc
Function to check whether a value is missing or not.
Space utilities, derived from NCAR software.
Definitions of constants and functions for working with missing values.
Distructor for triangles.
classe per la gestione del logging
Utilities for CHARACTER variables.