libsim Versione 7.2.1

◆ conxch_simc()

integer function conxch_simc ( double precision, dimension(:), intent(in) x,
double precision, dimension(:), intent(in) y,
integer, intent(in) i1,
integer, intent(in) i2,
integer, intent(in) i3,
integer, intent(in) i4 )
private

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.

Parametri
[in]yARRAYS CONTAINING THE COORDINATES OF THE DATA POINTS
[in]i1POINT NUMBERS OF FOUR POINTS P1,P2,P3, AND P4 THAT FORM A QUADRILATERAL WITH P3 AND P4 CONNECTED DIADONALLY.
[in]i2POINT NUMBERS OF FOUR POINTS P1,P2,P3, AND P4 THAT FORM A QUADRILATERAL WITH P3 AND P4 CONNECTED DIADONALLY.
[in]i3POINT NUMBERS OF FOUR POINTS P1,P2,P3, AND P4 THAT FORM A QUADRILATERAL WITH P3 AND P4 CONNECTED DIADONALLY.
[in]i4POINT NUMBERS OF FOUR POINTS P1,P2,P3, AND P4 THAT FORM A QUADRILATERAL WITH P3 AND P4 CONNECTED DIADONALLY.

Definizione alla linea 861 del file space_utilities.F90.

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

Generated with Doxygen.