libsim Versione 7.1.11
|
◆ conxch_simc()
THIS FUNCTION DETERMINES WHETHER OR NOT THE EXCHANGE OF TWO TRIANGLES IS NECESSARY ON THE BASIS OF MAX-MIN-ANGLE CRITERION BY C. L. LAWSON. THIS FUNCTION RETURNS A VALUE 1 (ONE) WHEN AN EXCHANGE IS NEEDED, AND 0 (ZERO) OTHERWISE.
Definizione alla linea 867 del file space_utilities.F90. 868! Copyright (C) 2011 ARPA-SIM <urpsim@smr.arpa.emr.it>
869! authors:
870! Davide Cesari <dcesari@arpa.emr.it>
871! Paolo Patruno <ppatruno@arpa.emr.it>
872
873! This program is free software; you can redistribute it and/or
874! modify it under the terms of the GNU General Public License as
875! published by the Free Software Foundation; either version 2 of
876! the License, or (at your option) any later version.
877
878! This program is distributed in the hope that it will be useful,
879! but WITHOUT ANY WARRANTY; without even the implied warranty of
880! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
881! GNU General Public License for more details.
882
883! You should have received a copy of the GNU General Public License
884! along with this program. If not, see <http://www.gnu.org/licenses/>.
885
886! This module contains contng derived from NCL - UNIX Version 6.0.0
887! $Id: contng.f,v 1.4 2008-07-27 00:16:57 haley Exp $
888!
889! Copyright (C) 2000
890! University Corporation for Atmospheric Research
891! All Rights Reserved
892!Redistribution and use in source and binary forms, with or without
893!modification, are permitted provided that the following conditions are
894!met:
895!
896!Neither the names of NCAR's Computational and Information Systems
897!Laboratory, the University Corporation for Atmospheric Research, nor
898!the names of its contributors may be used to endorse or promote
899!products derived from this Software without specific prior written
900!permission.
901!
902!Redistributions of source code must retain the above copyright
903!notice, this list of conditions, and the disclaimer below.
904!
905!Redistributions in binary form must reproduce the above copyright
906!notice, this list of conditions, and the disclaimer below in the
907!documentation and/or other materials provided with the distribution.
908!
909!THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
910!EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE WARRANTIES OF
911!MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
912!NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR COPYRIGHT
913!HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL,
914!EXEMPLARY, OR CONSEQUENTIAL DAMAGES OR OTHER LIABILITY, WHETHER IN AN
915!ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
916!CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS WITH THE
917!SOFTWARE.
918
919#include "config.h"
920
926
931implicit none
932
933type :: triangles
934 integer,pointer :: ipt(:) => null(), ipl(:) => null()
935 integer :: nt=imiss,nl=imiss
936end type triangles
937
938type :: xy
939 double precision :: x,y
940end type xy
941
945 MODULE PROCEDURE triangles_delete
946END INTERFACE
947
948INTERFACE triangles_compute
949 MODULE PROCEDURE triangles_compute_r, triangles_compute_d, triangles_compute_c
950END INTERFACE
951
952private
954
955
956contains
957
959function triangles_new(ndp) result(this)
960type(triangles) :: this
961integer,intent(in) :: ndp
962
963! those are done by type definition
964!this%nt=imiss
965!this%nl=imiss
966!nullify(this%ipt,this%ipl)
967
969 allocate(this%ipt(6*ndp-15), this%ipl(6*ndp))
970else
971 this%nt=0
972 this%nl=0
973end if
974return
975end function triangles_new
976
977
979subroutine triangles_delete(this)
980type(triangles) :: this
981
982if (associated(this%ipt)) deallocate(this%ipt)
983if (associated(this%ipl)) deallocate(this%ipl)
984
985this%nt=imiss
986this%nl=imiss
987
988end subroutine triangles_delete
989
990
991integer function triangles_compute_r (XD,YD,tri)
992real,intent(in) :: XD(:)
993real,intent(in) :: YD(:)
994type (triangles),intent(inout) :: tri
995type (xy) :: co(size(xd))
996
997if (tri%nt /= 0) then
998 co%x=dble(xd)
999 co%y=dble(yd)
1000 triangles_compute_r = contng_simc(co,tri%NT,tri%IPT,tri%NL,tri%IPL)
1001end if
1002end function triangles_compute_r
1003
1004integer function triangles_compute_d (XD,YD,tri)
1005double precision,intent(in) :: XD(:)
1006double precision,intent(in) :: YD(:)
1007type (triangles),intent(inout) :: tri
1008type (xy) :: co(size(xd))
1009
1010if (tri%nt /= 0) then
1011 co%x=xd
1012 co%y=yd
1013 triangles_compute_d = contng_simc(co,tri%NT,tri%IPT,tri%NL,tri%IPL)
1014end if
1015end function triangles_compute_d
1016
1017integer function triangles_compute_c (co,tri)
1018type (xy),intent(in) :: co(:)
1019type (triangles),intent(inout) :: tri
1020
1021if (tri%nt /= 0) then
1022 triangles_compute_c = contng_simc(co,tri%NT,tri%IPT,tri%NL,tri%IPL)
1023end if
1024end function triangles_compute_c
1025
1026
1040integer function contng_simc (co,NT,IPT,NL,IPL)
1041
1042type (xy), intent(in) :: co(:)
1043integer, intent(out) :: NT
1044integer, intent(out) :: NL
1049integer,intent(out) :: IPT(:)
1055integer,intent(out) :: IPL(:)
1056
1057!!$C THE internal PARAMETERS ARE
1058!!$C IWL = INTEGER ARRAY OF DIMENSION 18*NDP USED
1059!!$C INTERNALLY AS A WORK AREA,
1060!!$C IWP = INTEGER ARRAY OF DIMENSION NDP USED
1061!!$C INTERNALLY AS A WORK AREA,
1062!!$C WK = ARRAY OF DIMENSION NDP USED INTERNALLY AS A
1063!!$C WORK AREA.
1064integer :: IWL(18*size(co)), IWP(size(co))
1065
1066double precision :: WK(size(co))
1067integer :: ITF(2), NREP=100
1068real :: RATIO=1.0e-6
1069double precision :: AR,ARMN,ARMX,DSQ12,DSQI,DSQMN,DSQMX,DX,DX21,DXMN,DXMX,DY,DY21,DYMN,DYMX
1070double precision :: X1,Y1
1071integer :: ilf,IP,IP1,IP1P1,IP2,IP3,IPL1,IPL2,IPLJ1,IPLJ2,IPMN1,IPMN2,IPT1,IPT2,IPT3
1072INTEGER :: IPTI,IPTI1,IPTI2,IREP,IT,IT1T3,IT2T3,ITS,ITT3
1073INTEGER :: ILFT2,ITT3R,JLT3,JP,JP1,JP2,JP2T3,JP3T3,JPC,JPMN,JPMX,JWL,JWL1
1074INTEGER :: JWL1MN,NDP,NDPM1,NL0,NLF,NLFC,NLFT2,NLN,NLNT3,NLT3,NSH,NSHT3,NT0,NTF
1075INTEGER :: NTT3,NTT3P3
1076logical :: err
1077
1078integer ::i,mloc(size(co)-1),mlocall(1),mlocv(1)
1079type(xy) :: dmp
1080
1081 !
1082 ! PRELIMINARY PROCESSING
1083 !
1084
1085nt = imiss
1086nl = imiss
1087
1088contng_simc=0
1089ndp=size(co)
1090ndpm1 = ndp-1
1091 !
1092 ! DETERMINES THE CLOSEST PAIR OF DATA POINTS AND THEIR MIDPOINT.
1093 !
1094
1095call l4f_log(l4f_debug,"start triangulation")
1096
1097do i=1,size(co)-1
1098 mlocv=minloc(vdsqf(co(i),co(i+1:)))+i
1099 mloc(i)=mlocv(1)
1100end do
1101
1102mlocall=minloc((/(vdsqf(co(i),co(mloc(i))),i=1,size(mloc))/))
1103
1104dsqmn = vdsqf(co(mlocall(1)),co(mloc(mlocall(1))))
1105ipmn1 = mlocall(1)
1106ipmn2 = mloc(mlocall(1))
1107
1108call l4f_log(l4f_debug,"end triangulation closest pair")
1109!!$print *, DSQMN, IPMN1, IPMN2
1110
1111IF (dsqmn == 0.) then
1112 !
1113 ! ERROR, IDENTICAL INPUT DATA POINTS
1114 !
1115 call l4f_log(l4f_error,"CONTNG-IDENTICAL INPUT DATA POINTS FOUND")
1116 contng_simc=1
1117 RETURN
1118end IF
1119
1120!!$call l4f_log(L4F_DEBUG,"start your")
1121!!$
1122!!$DSQMN = DSQF(XD(1),YD(1),XD(2),YD(2))
1123!!$IPMN1 = 1
1124!!$IPMN2 = 2
1125!!$DO IP1=1,NDPM1
1126!!$ X1 = XD(IP1)
1127!!$ Y1 = YD(IP1)
1128!!$ IP1P1 = IP1+1
1129!!$ DO IP2=IP1P1,NDP
1130!!$ DSQI = DSQF(X1,Y1,XD(IP2),YD(IP2))
1131!!$
1132!!$ IF (DSQI == 0.) then
1133!!$ !
1134!!$ ! ERROR, IDENTICAL INPUT DATA POINTS
1135!!$ !
1136!!$ call l4f_log(L4F_ERROR,"CONTNG-IDENTICAL INPUT DATA POINTS FOUND AT "//t2c(ip1)//" AND"//t2c(ip2))
1137!!$ CONTNG_simc=1
1138!!$ RETURN
1139!!$ end IF
1140!!$
1141!!$ IF (DSQI .GE. DSQMN) CYCLE
1142!!$ DSQMN = DSQI
1143!!$ IPMN1 = IP1
1144!!$ IPMN2 = IP2
1145!!$ end DO
1146!!$end DO
1147!!$
1148!!$call l4f_log(L4F_DEBUG,"end your")
1149!!$print *, DSQMN, IPMN1, IPMN2
1150
1151dsq12 = dsqmn
1152dmp%x = (co(ipmn1)%x+co(ipmn2)%x)/2.0
1153dmp%y = (co(ipmn1)%y+co(ipmn2)%y)/2.0
1154 !
1155 ! SORTS THE OTHER (NDP-2) DATA POINTS IN ASCENDING ORDER OF
1156 ! DISTANCE FROM THE MIDPOINT AND STORES THE SORTED DATA POINT
1157 ! NUMBERS IN THE IWP ARRAY.
1158 !
1159jp1 = 2
1160DO ip1=1,ndp
1161 IF (ip1.EQ.ipmn1 .OR. ip1.EQ.ipmn2) cycle
1162 jp1 = jp1+1
1163 iwp(jp1) = ip1
1164 wk(jp1) = vdsqf(dmp,co(ip1))
1165end DO
1166DO jp1=3,ndpm1
1167 dsqmn = wk(jp1)
1168 jpmn = jp1
1169 DO jp2=jp1,ndp
1170 IF (wk(jp2) .GE. dsqmn) cycle
1171 dsqmn = wk(jp2)
1172 jpmn = jp2
1173 end DO
1174 its = iwp(jp1)
1175 iwp(jp1) = iwp(jpmn)
1176 iwp(jpmn) = its
1177 wk(jpmn) = wk(jp1)
1178end DO
1179
1180call l4f_log(l4f_debug,"end triangulation sort")
1181
1182 !
1183 ! IF NECESSARY, MODIFIES THE ORDERING IN SUCH A WAY THAT THE
1184 ! FIRST THREE DATA POINTS ARE NOT COLLINEAR.
1185 !
1186ar = dsq12*ratio
1187x1 = co(ipmn1)%x
1188y1 = co(ipmn1)%y
1189dx21 = co(ipmn2)%x-x1
1190dy21 = co(ipmn2)%y-y1
1191
1192err=.true.
1193DO jp=3,ndp
1194 ip = iwp(jp)
1195 IF (abs((co(ip)%y-y1)*dx21-(co(ip)%x-x1)*dy21) .GT. ar) then
1196 err=.false.
1197 exit
1198 end IF
1199end DO
1200if (err) then
1201 call l4f_log(l4f_debug,"CONTNG - ALL COLLINEAR DATA POINTS")
1202 contng_simc=2
1203 return
1204end if
1205IF (jp /= 3) then
1206 jpmx = jp
1207 jp = jpmx+1
1208 DO jpc=4,jpmx
1209 jp = jp-1
1210 iwp(jp) = iwp(jp-1)
1211 end DO
1212 iwp(3) = ip
1213end IF
1214call l4f_log(l4f_debug,"end triangulation collinear")
1215
1216 !
1217 ! FORMS THE FIRST TRIANGLE. STORES POINT NUMBERS OF THE VER-
1218 ! TEXES OF THE TRIANGLE IN THE IPT ARRAY, AND STORES POINT NUM-
1219 ! BERS OF THE BORDER LINE SEGMENTS AND THE TRIANGLE NUMBER IN
1220 ! THE IPL ARRAY.
1221 !
1222ip1 = ipmn1
1223ip2 = ipmn2
1224ip3 = iwp(3)
1225IF (side(co(ip1),co(ip2),co(ip3)) < 10.0) then
1226 ip1 = ipmn2
1227 ip2 = ipmn1
1228end IF
1229
1230nt0 = 1
1231ntt3 = 3
1232ipt(1) = ip1
1233ipt(2) = ip2
1234ipt(3) = ip3
1235nl0 = 3
1236nlt3 = 9
1237ipl(1) = ip1
1238ipl(2) = ip2
1239ipl(3) = 1
1240ipl(4) = ip2
1241ipl(5) = ip3
1242ipl(6) = 1
1243ipl(7) = ip3
1244ipl(8) = ip1
1245ipl(9) = 1
1246
1247call l4f_log(l4f_debug,"end triangulation first triangle")
1248
1249 !
1250 ! ADDS THE REMAINING (NDP-3) DATA POINTS, ONE BY ONE.
1251 !
1252l400 : DO jp1=4,ndp
1253 ip1 = iwp(jp1)
1254 x1 = co(ip1)%x
1255 y1 = co(ip1)%y
1256 !
1257 ! - DETERMINES THE VISIBLE BORDER LINE SEGMENTS.
1258 !
1259 ip2 = ipl(1)
1260 jpmn = 1
1261 dxmn = co(ip2)%x-x1
1262 dymn = co(ip2)%y-y1
1263 dsqmn = dxmn**2+dymn**2
1264 armn = dsqmn*ratio
1265 jpmx = 1
1266 dxmx = dxmn
1267 dymx = dymn
1268 dsqmx = dsqmn
1269 armx = armn
1270 DO jp2=2,nl0
1271 ip2 = ipl(3*jp2-2)
1272 dx = co(ip2)%x-x1
1273 dy = co(ip2)%y-y1
1274 ar = dy*dxmn-dx*dymn
1275 IF (ar <= armn) then
1276 dsqi = dx**2+dy**2
1277 IF (ar.GE.(-armn) .AND. dsqi.GE.dsqmn) GO TO 230
1278 jpmn = jp2
1279 dxmn = dx
1280 dymn = dy
1281 dsqmn = dsqi
1282 armn = dsqmn*ratio
1283 end IF
1284230 ar = dy*dxmx-dx*dymx
1285 IF (ar .LT. (-armx)) cycle
1286 dsqi = dx**2+dy**2
1287 IF (ar.LE.armx .AND. dsqi.GE.dsqmx) cycle
1288 jpmx = jp2
1289 dxmx = dx
1290 dymx = dy
1291 dsqmx = dsqi
1292 armx = dsqmx*ratio
1293 end DO
1294 IF (jpmx .LT. jpmn) jpmx = jpmx+nl0
1295 nsh = jpmn-1
1296 IF (nsh > 0) then
1297 !
1298 ! - SHIFTS (ROTATES) THE IPL ARRAY TO HAVE THE INVISIBLE BORDER
1299 ! - LINE SEGMENTS CONTAINED IN THE FIRST PART OF THE IPL ARRAY.
1300 !
1301 nsht3 = nsh*3
1302 DO jp2t3=3,nsht3,3
1303 jp3t3 = jp2t3+nlt3
1304 ipl(jp3t3-2) = ipl(jp2t3-2)
1305 ipl(jp3t3-1) = ipl(jp2t3-1)
1306 ipl(jp3t3) = ipl(jp2t3)
1307 end DO
1308 DO jp2t3=3,nlt3,3
1309 jp3t3 = jp2t3+nsht3
1310 ipl(jp2t3-2) = ipl(jp3t3-2)
1311 ipl(jp2t3-1) = ipl(jp3t3-1)
1312 ipl(jp2t3) = ipl(jp3t3)
1313 end DO
1314 jpmx = jpmx-nsh
1315 !
1316 ! - ADDS TRIANGLES TO THE IPT ARRAY, UPDATES BORDER LINE
1317 ! - SEGMENTS IN THE IPL ARRAY, AND SETS FLAGS FOR THE BORDER
1318 ! - LINE SEGMENTS TO BE REEXAMINED IN THE IWL ARRAY.
1319 !
1320 end IF
1321
1322 jwl = 0
1323 l310 : DO jp2=jpmx,nl0
1324 jp2t3 = jp2*3
1325 ipl1 = ipl(jp2t3-2)
1326 ipl2 = ipl(jp2t3-1)
1327 it = ipl(jp2t3)
1328 !
1329 ! - - ADDS A TRIANGLE TO THE IPT ARRAY.
1330 !
1331 nt0 = nt0+1
1332 ntt3 = ntt3+3
1333 ipt(ntt3-2) = ipl2
1334 ipt(ntt3-1) = ipl1
1335 ipt(ntt3) = ip1
1336 !
1337 ! - - UPDATES BORDER LINE SEGMENTS IN THE IPL ARRAY.
1338 !
1339 IF (jp2 == jpmx) then
1340 ipl(jp2t3-1) = ip1
1341 ipl(jp2t3) = nt0
1342 end IF
1343 IF (jp2 == nl0) then
1344 nln = jpmx+1
1345 nlnt3 = nln*3
1346 ipl(nlnt3-2) = ip1
1347 ipl(nlnt3-1) = ipl(1)
1348 ipl(nlnt3) = nt0
1349 end IF
1350 !
1351 ! - - DETERMINES THE VERTEX THAT DOES NOT LIE ON THE BORDER
1352 ! - - LINE SEGMENTS.
1353 !
1354 itt3 = it*3
1355 ipti = ipt(itt3-2)
1356 IF (ipti.NE.ipl1 .AND. ipti.NE.ipl2) GO TO 300
1357 ipti = ipt(itt3-1)
1358 IF (ipti.NE.ipl1 .AND. ipti.NE.ipl2) GO TO 300
1359 ipti = ipt(itt3)
1360 !
1361 ! - - CHECKS IF THE EXCHANGE IS NECESSARY.
1362 !
1363300 IF (conxch_simc(co%X,co%Y,ip1,ipti,ipl1,ipl2) .EQ. 0) cycle l310
1364 !
1365 ! - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
1366 !
1367 ipt(itt3-2) = ipti
1368 ipt(itt3-1) = ipl1
1369 ipt(itt3) = ip1
1370 ipt(ntt3-1) = ipti
1371 IF (jp2 .EQ. jpmx) ipl(jp2t3) = it
1372 IF (jp2.EQ.nl0 .AND. ipl(3).EQ.it) ipl(3) = nt0
1373 !
1374 ! - - SETS FLAGS IN THE IWL ARRAY.
1375 !
1376 jwl = jwl+4
1377 iwl(jwl-3) = ipl1
1378 iwl(jwl-2) = ipti
1379 iwl(jwl-1) = ipti
1380 iwl(jwl) = ipl2
1381 end DO l310
1382 nl0 = nln
1383 nlt3 = nlnt3
1384 nlf = jwl/2
1385 IF (nlf .EQ. 0) cycle l400
1386 !
1387 ! - IMPROVES TRIANGULATION.
1388 !
1389 ntt3p3 = ntt3+3
1390 DO irep=1,nrep
1391 l370 : DO ilf=1,nlf
1392 ilft2 = ilf*2
1393 ipl1 = iwl(ilft2-1)
1394 ipl2 = iwl(ilft2)
1395 !
1396 ! - - LOCATES IN THE IPT ARRAY TWO TRIANGLES ON BOTH SIDES OF
1397 ! - - THE FLAGGED LINE SEGMENT.
1398 !
1399 ntf = 0
1400 DO itt3r=3,ntt3,3
1401 itt3 = ntt3p3-itt3r
1402 ipt1 = ipt(itt3-2)
1403 ipt2 = ipt(itt3-1)
1404 ipt3 = ipt(itt3)
1405 IF (ipl1.NE.ipt1 .AND. ipl1.NE.ipt2 .AND. ipl1.NE.ipt3) cycle
1406 IF (ipl2.NE.ipt1 .AND. ipl2.NE.ipt2 .AND. ipl2.NE.ipt3) cycle
1407 ntf = ntf+1
1408 itf(ntf) = itt3/3
1409 IF (ntf .EQ. 2) GO TO 330
1410 end DO
1411 IF (ntf .LT. 2) cycle l370
1412 !
1413 ! - - DETERMINES THE VERTEXES OF THE TRIANGLES THAT DO NOT LIE
1414 ! - - ON THE LINE SEGMENT.
1415 !
1416330 it1t3 = itf(1)*3
1417 ipti1 = ipt(it1t3-2)
1418 IF (ipti1.NE.ipl1 .AND. ipti1.NE.ipl2) GO TO 340
1419 ipti1 = ipt(it1t3-1)
1420 IF (ipti1.NE.ipl1 .AND. ipti1.NE.ipl2) GO TO 340
1421 ipti1 = ipt(it1t3)
1422340 it2t3 = itf(2)*3
1423 ipti2 = ipt(it2t3-2)
1424 IF (ipti2.NE.ipl1 .AND. ipti2.NE.ipl2) GO TO 350
1425 ipti2 = ipt(it2t3-1)
1426 IF (ipti2.NE.ipl1 .AND. ipti2.NE.ipl2) GO TO 350
1427 ipti2 = ipt(it2t3)
1428 !
1429 ! - - CHECKS IF THE EXCHANGE IS NECESSARY.
1430 !
1431350 IF (conxch_simc(co%X,co%Y,ipti1,ipti2,ipl1,ipl2) .EQ. 0) cycle l370
1432 !
1433 ! - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
1434 !
1435 ipt(it1t3-2) = ipti1
1436 ipt(it1t3-1) = ipti2
1437 ipt(it1t3) = ipl1
1438 ipt(it2t3-2) = ipti2
1439 ipt(it2t3-1) = ipti1
1440 ipt(it2t3) = ipl2
1441 !
1442 ! - - SETS NEW FLAGS.
1443 !
1444 jwl = jwl+8
1445 iwl(jwl-7) = ipl1
1446 iwl(jwl-6) = ipti1
1447 iwl(jwl-5) = ipti1
1448 iwl(jwl-4) = ipl2
1449 iwl(jwl-3) = ipl2
1450 iwl(jwl-2) = ipti2
1451 iwl(jwl-1) = ipti2
1452 iwl(jwl) = ipl1
1453 DO jlt3=3,nlt3,3
1454 iplj1 = ipl(jlt3-2)
1455 iplj2 = ipl(jlt3-1)
1456 IF ((iplj1.EQ.ipl1 .AND. iplj2.EQ.ipti2) .OR. &
1457 (iplj2.EQ.ipl1 .AND. iplj1.EQ.ipti2)) ipl(jlt3) = itf(1)
1458 IF ((iplj1.EQ.ipl2 .AND. iplj2.EQ.ipti1) .OR. &
1459 (iplj2.EQ.ipl2 .AND. iplj1.EQ.ipti1)) ipl(jlt3) = itf(2)
1460 end DO
1461 end DO l370
1462 nlfc = nlf
1463 nlf = jwl/2
1464 IF (nlf .EQ. nlfc) cycle l400
1465 !
1466 ! - - RESETS THE IWL ARRAY FOR THE NEXT ROUND.
1467 !
1468 jwl = 0
1469 jwl1mn = (nlfc+1)*2
1470 nlft2 = nlf*2
1471 DO jwl1=jwl1mn,nlft2,2
1472 jwl = jwl+2
1473 iwl(jwl-1) = iwl(jwl1-1)
1474 iwl(jwl) = iwl(jwl1)
1475 end DO
1476 nlf = jwl/2
1477 end DO
1478end DO l400
1479
1480call l4f_log(l4f_debug,"end triangulation appending")
1481
1482 !
1483 ! REARRANGE THE IPT ARRAY SO THAT THE VERTEXES OF EACH TRIANGLE
1484 ! ARE LISTED COUNTER-CLOCKWISE.
1485 !
1486DO itt3=3,ntt3,3
1487 ip1 = ipt(itt3-2)
1488 ip2 = ipt(itt3-1)
1489 ip3 = ipt(itt3)
1490 IF (side(co(ip1),co(ip2),co(ip3)) .GE. 10.0) cycle
1491 ipt(itt3-2) = ip2
1492 ipt(itt3-1) = ip1
1493end DO
1494call l4f_log(l4f_debug,"end triangulation rearranging")
1495
1496nt = nt0
1497nl = nl0
1498
1499call l4f_log(l4f_debug,"end triangulation")
1500
1501RETURN
1502
1503contains
1504
1505!!$double precision function DSQF(U1,V1,U2,V2)
1506!!$double precision,intent(in) :: U1,V1,U2,V2
1507!!$
1508!!$DSQF = (U2-U1)**2+(V2-V1)**2
1509!!$end function DSQF
1510
1511
1512elemental double precision function vdsqf(co1,co2)
1513type(xy),intent(in) :: co1,co2
1514
1515vdsqf = (co2%x-co1%x)**2+(co2%y-co1%y)**2
1516if (vdsqf == 0.d0) vdsqf = huge(vdsqf)
1517end function vdsqf
1518
1519
1520!!$double precision function SIDE(U1,V1,U2,V2,U3,V3)
1521!!$double precision,intent(in):: U1,V1,U2,V2,U3,V3
1522!!$
1523!!$SIDE = (V3-V1)*(U2-U1)-(U3-U1)*(V2-V1)
1524!!$end function SIDE
1525
1526double precision function side(co1,co2,co3)
1527type(xy),intent(in):: co1,co2,co3
1528
1529side = (co3%y-co1%y)*(co2%x-co1%x)-(co3%x-co1%x)*(co2%y-co1%y)
1530end function side
1531
1532end function contng_simc
1533
1534
1540INTEGER FUNCTION conxch_simc (X,Y,I1,I2,I3,I4)
1541
1542double precision,intent(in) :: X(:), Y(:)
1545integer,intent(in) :: I1,I2,I3,I4
1546
1547double precision :: X0(4), Y0(4)
1548double precision :: C2SQ,C1SQ,A3SQ,B2SQ,B3SQ,A1SQ,A4SQ,B1SQ,B4SQ,A2SQ,C4SQ,C3SQ
1549integer :: IDX
1550double precision :: S1SQ,S2SQ,S3SQ,S4SQ,U1,U2,U3,U4
1551
1552equivalence(c2sq,c1sq),(a3sq,b2sq),(b3sq,a1sq),(a4sq,b1sq), &
1553 (b4sq,a2sq),(c4sq,c3sq)
1554
1555x0(1) = x(i1)
1556y0(1) = y(i1)
1557x0(2) = x(i2)
1558y0(2) = y(i2)
1559x0(3) = x(i3)
1560y0(3) = y(i3)
1561x0(4) = x(i4)
1562y0(4) = y(i4)
1563idx = 0
1564u3 = (y0(2)-y0(3))*(x0(1)-x0(3))-(x0(2)-x0(3))*(y0(1)-y0(3))
1565u4 = (y0(1)-y0(4))*(x0(2)-x0(4))-(x0(1)-x0(4))*(y0(2)-y0(4))
1566IF (u3*u4 > 0.0) then
1567 u1 = (y0(3)-y0(1))*(x0(4)-x0(1))-(x0(3)-x0(1))*(y0(4)-y0(1))
1568 u2 = (y0(4)-y0(2))*(x0(3)-x0(2))-(x0(4)-x0(2))*(y0(3)-y0(2))
1569 a1sq = (x0(1)-x0(3))**2+(y0(1)-y0(3))**2
1570 b1sq = (x0(4)-x0(1))**2+(y0(4)-y0(1))**2
1571 c1sq = (x0(3)-x0(4))**2+(y0(3)-y0(4))**2
1572 a2sq = (x0(2)-x0(4))**2+(y0(2)-y0(4))**2
1573 b2sq = (x0(3)-x0(2))**2+(y0(3)-y0(2))**2
1574 c3sq = (x0(2)-x0(1))**2+(y0(2)-y0(1))**2
1575 s1sq = u1*u1/(c1sq*max(a1sq,b1sq))
1576 s2sq = u2*u2/(c2sq*max(a2sq,b2sq))
1577 s3sq = u3*u3/(c3sq*max(a3sq,b3sq))
1578 s4sq = u4*u4/(c4sq*max(a4sq,b4sq))
1579 IF (min(s1sq,s2sq) < min(s3sq,s4sq)) idx = 1
1580end IF
1581conxch_simc = idx
1582RETURN
1583
1584END FUNCTION conxch_simc
1585
Function to check whether a value is missing or not. Definition: missing_values.f90:72 Definitions of constants and functions for working with missing values. Definition: missing_values.f90:50 |