libsim Versione 7.2.1
space_utilities.F90
1! Copyright (C) 2011 ARPA-SIM <urpsim@smr.arpa.emr.it>
2! authors:
3! Davide Cesari <dcesari@arpa.emr.it>
4! Paolo Patruno <ppatruno@arpa.emr.it>
5
6! This program is free software; you can redistribute it and/or
7! modify it under the terms of the GNU General Public License as
8! published by the Free Software Foundation; either version 2 of
9! the License, or (at your option) any later version.
10
11! This program is distributed in the hope that it will be useful,
12! but WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14! GNU General Public License for more details.
15
16! You should have received a copy of the GNU General Public License
17! along with this program. If not, see <http://www.gnu.org/licenses/>.
18
19! This module contains contng derived from NCL - UNIX Version 6.0.0
20! $Id: contng.f,v 1.4 2008-07-27 00:16:57 haley Exp $
21!
22! Copyright (C) 2000
23! University Corporation for Atmospheric Research
24! All Rights Reserved
25!Redistribution and use in source and binary forms, with or without
26!modification, are permitted provided that the following conditions are
27!met:
28!
29!Neither the names of NCAR's Computational and Information Systems
30!Laboratory, the University Corporation for Atmospheric Research, nor
31!the names of its contributors may be used to endorse or promote
32!products derived from this Software without specific prior written
33!permission.
34!
35!Redistributions of source code must retain the above copyright
36!notice, this list of conditions, and the disclaimer below.
37!
38!Redistributions in binary form must reproduce the above copyright
39!notice, this list of conditions, and the disclaimer below in the
40!documentation and/or other materials provided with the distribution.
41!
42!THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
43!EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE WARRANTIES OF
44!MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
45!NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR COPYRIGHT
46!HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL,
47!EXEMPLARY, OR CONSEQUENTIAL DAMAGES OR OTHER LIABILITY, WHETHER IN AN
48!ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
49!CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS WITH THE
50!SOFTWARE.
51
52#include "config.h"
53
59
60module space_utilities
64implicit none
65
66type :: triangles
67 integer,pointer :: ipt(:) => null(), ipl(:) => null()
68 integer :: nt=imiss,nl=imiss
69end type triangles
70
71type :: xy
72 double precision :: x,y
73end type xy
74
77INTERFACE delete
78 MODULE PROCEDURE triangles_delete
79END INTERFACE
80
81INTERFACE triangles_compute
82 MODULE PROCEDURE triangles_compute_r, triangles_compute_d, triangles_compute_c
83END INTERFACE
84
85private
86public triangles, triangles_new, delete, triangles_compute, xy
87
88
89contains
90
92function triangles_new(ndp) result(this)
93type(triangles) :: this
94integer,intent(in) :: ndp
95
96! those are done by type definition
97!this%nt=imiss
98!this%nl=imiss
99!nullify(this%ipt,this%ipl)
100
101if (c_e(ndp) .and. ndp >= 3)then
102 allocate(this%ipt(6*ndp-15), this%ipl(6*ndp))
103else
104 this%nt=0
105 this%nl=0
106end if
107return
108end function triangles_new
109
110
112subroutine triangles_delete(this)
113type(triangles) :: this
114
115if (associated(this%ipt)) deallocate(this%ipt)
116if (associated(this%ipl)) deallocate(this%ipl)
117
118this%nt=imiss
119this%nl=imiss
120
121end subroutine triangles_delete
122
123
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))
129
130if (tri%nt /= 0) then
131 co%x=dble(xd)
132 co%y=dble(yd)
133 triangles_compute_r = contng_simc(co,tri%NT,tri%IPT,tri%NL,tri%IPL)
134end if
135end function triangles_compute_r
136
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))
142
143if (tri%nt /= 0) then
144 co%x=xd
145 co%y=yd
146 triangles_compute_d = contng_simc(co,tri%NT,tri%IPT,tri%NL,tri%IPL)
147end if
148end function triangles_compute_d
149
150integer function triangles_compute_c (co,tri)
151type (xy),intent(in) :: co(:)
152type (triangles),intent(inout) :: tri
153
154if (tri%nt /= 0) then
155 triangles_compute_c = contng_simc(co,tri%NT,tri%IPT,tri%NL,tri%IPL)
156end if
157end function triangles_compute_c
158
159
173integer function contng_simc (co,NT,IPT,NL,IPL)
174
175type (xy), intent(in) :: co(:)
176integer, intent(out) :: NT
177integer, intent(out) :: NL
182integer,intent(out) :: IPT(:)
188integer,intent(out) :: IPL(:)
189
190!!$C THE internal PARAMETERS ARE
191!!$C IWL = INTEGER ARRAY OF DIMENSION 18*NDP USED
192!!$C INTERNALLY AS A WORK AREA,
193!!$C IWP = INTEGER ARRAY OF DIMENSION NDP USED
194!!$C INTERNALLY AS A WORK AREA,
195!!$C WK = ARRAY OF DIMENSION NDP USED INTERNALLY AS A
196!!$C WORK AREA.
197integer :: IWL(18*size(co)), IWP(size(co))
198
199double precision :: WK(size(co))
200integer :: ITF(2), NREP=100
201real :: RATIO=1.0e-6
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
209logical :: err
210
211integer ::i,mloc(size(co)-1),mlocall(1),mlocv(1)
212type(xy) :: dmp
213
214 !
215 ! PRELIMINARY PROCESSING
216 !
217
218nt = imiss
219nl = imiss
220
221contng_simc=0
222ndp=size(co)
223ndpm1 = ndp-1
224 !
225 ! DETERMINES THE CLOSEST PAIR OF DATA POINTS AND THEIR MIDPOINT.
226 !
227
228call l4f_log(l4f_debug,"start triangulation")
229
230do i=1,size(co)-1
231 mlocv=minloc(vdsqf(co(i),co(i+1:)))+i
232 mloc(i)=mlocv(1)
233end do
234
235mlocall=minloc((/(vdsqf(co(i),co(mloc(i))),i=1,size(mloc))/))
236
237dsqmn = vdsqf(co(mlocall(1)),co(mloc(mlocall(1))))
238ipmn1 = mlocall(1)
239ipmn2 = mloc(mlocall(1))
240
241call l4f_log(l4f_debug,"end triangulation closest pair")
242!!$print *, DSQMN, IPMN1, IPMN2
243
244IF (dsqmn == 0.) then
245 !
246 ! ERROR, IDENTICAL INPUT DATA POINTS
247 !
248 call l4f_log(l4f_error,"CONTNG-IDENTICAL INPUT DATA POINTS FOUND")
249 contng_simc=1
250 RETURN
251end IF
252
253!!$call l4f_log(L4F_DEBUG,"start your")
254!!$
255!!$DSQMN = DSQF(XD(1),YD(1),XD(2),YD(2))
256!!$IPMN1 = 1
257!!$IPMN2 = 2
258!!$DO IP1=1,NDPM1
259!!$ X1 = XD(IP1)
260!!$ Y1 = YD(IP1)
261!!$ IP1P1 = IP1+1
262!!$ DO IP2=IP1P1,NDP
263!!$ DSQI = DSQF(X1,Y1,XD(IP2),YD(IP2))
264!!$
265!!$ IF (DSQI == 0.) then
266!!$ !
267!!$ ! ERROR, IDENTICAL INPUT DATA POINTS
268!!$ !
269!!$ call l4f_log(L4F_ERROR,"CONTNG-IDENTICAL INPUT DATA POINTS FOUND AT "//t2c(ip1)//" AND"//t2c(ip2))
270!!$ CONTNG_simc=1
271!!$ RETURN
272!!$ end IF
273!!$
274!!$ IF (DSQI .GE. DSQMN) CYCLE
275!!$ DSQMN = DSQI
276!!$ IPMN1 = IP1
277!!$ IPMN2 = IP2
278!!$ end DO
279!!$end DO
280!!$
281!!$call l4f_log(L4F_DEBUG,"end your")
282!!$print *, DSQMN, IPMN1, IPMN2
283
284dsq12 = dsqmn
285dmp%x = (co(ipmn1)%x+co(ipmn2)%x)/2.0
286dmp%y = (co(ipmn1)%y+co(ipmn2)%y)/2.0
287 !
288 ! SORTS THE OTHER (NDP-2) DATA POINTS IN ASCENDING ORDER OF
289 ! DISTANCE FROM THE MIDPOINT AND STORES THE SORTED DATA POINT
290 ! NUMBERS IN THE IWP ARRAY.
291 !
292jp1 = 2
293DO ip1=1,ndp
294 IF (ip1.EQ.ipmn1 .OR. ip1.EQ.ipmn2) cycle
295 jp1 = jp1+1
296 iwp(jp1) = ip1
297 wk(jp1) = vdsqf(dmp,co(ip1))
298end DO
299DO jp1=3,ndpm1
300 dsqmn = wk(jp1)
301 jpmn = jp1
302 DO jp2=jp1,ndp
303 IF (wk(jp2) .GE. dsqmn) cycle
304 dsqmn = wk(jp2)
305 jpmn = jp2
306 end DO
307 its = iwp(jp1)
308 iwp(jp1) = iwp(jpmn)
309 iwp(jpmn) = its
310 wk(jpmn) = wk(jp1)
311end DO
313call l4f_log(l4f_debug,"end triangulation sort")
314
315 !
316 ! IF NECESSARY, MODIFIES THE ORDERING IN SUCH A WAY THAT THE
317 ! FIRST THREE DATA POINTS ARE NOT COLLINEAR.
318 !
319ar = dsq12*ratio
320x1 = co(ipmn1)%x
321y1 = co(ipmn1)%y
322dx21 = co(ipmn2)%x-x1
323dy21 = co(ipmn2)%y-y1
324
325err=.true.
326DO jp=3,ndp
327 ip = iwp(jp)
328 IF (abs((co(ip)%y-y1)*dx21-(co(ip)%x-x1)*dy21) .GT. ar) then
329 err=.false.
330 exit
331 end IF
332end DO
333if (err) then
334 call l4f_log(l4f_debug,"CONTNG - ALL COLLINEAR DATA POINTS")
335 contng_simc=2
336 return
337end if
338IF (jp /= 3) then
339 jpmx = jp
340 jp = jpmx+1
341 DO jpc=4,jpmx
342 jp = jp-1
343 iwp(jp) = iwp(jp-1)
344 end DO
345 iwp(3) = ip
346end IF
347call l4f_log(l4f_debug,"end triangulation collinear")
348
349 !
350 ! FORMS THE FIRST TRIANGLE. STORES POINT NUMBERS OF THE VER-
351 ! TEXES OF THE TRIANGLE IN THE IPT ARRAY, AND STORES POINT NUM-
352 ! BERS OF THE BORDER LINE SEGMENTS AND THE TRIANGLE NUMBER IN
353 ! THE IPL ARRAY.
354 !
355ip1 = ipmn1
356ip2 = ipmn2
357ip3 = iwp(3)
358IF (side(co(ip1),co(ip2),co(ip3)) < 10.0) then
359 ip1 = ipmn2
360 ip2 = ipmn1
361end IF
362
363nt0 = 1
364ntt3 = 3
365ipt(1) = ip1
366ipt(2) = ip2
367ipt(3) = ip3
368nl0 = 3
369nlt3 = 9
370ipl(1) = ip1
371ipl(2) = ip2
372ipl(3) = 1
373ipl(4) = ip2
374ipl(5) = ip3
375ipl(6) = 1
376ipl(7) = ip3
377ipl(8) = ip1
378ipl(9) = 1
379
380call l4f_log(l4f_debug,"end triangulation first triangle")
381
382 !
383 ! ADDS THE REMAINING (NDP-3) DATA POINTS, ONE BY ONE.
384 !
385l400 : DO jp1=4,ndp
386 ip1 = iwp(jp1)
387 x1 = co(ip1)%x
388 y1 = co(ip1)%y
389 !
390 ! - DETERMINES THE VISIBLE BORDER LINE SEGMENTS.
391 !
392 ip2 = ipl(1)
393 jpmn = 1
394 dxmn = co(ip2)%x-x1
395 dymn = co(ip2)%y-y1
396 dsqmn = dxmn**2+dymn**2
397 armn = dsqmn*ratio
398 jpmx = 1
399 dxmx = dxmn
400 dymx = dymn
401 dsqmx = dsqmn
402 armx = armn
403 DO jp2=2,nl0
404 ip2 = ipl(3*jp2-2)
405 dx = co(ip2)%x-x1
406 dy = co(ip2)%y-y1
407 ar = dy*dxmn-dx*dymn
408 IF (ar <= armn) then
409 dsqi = dx**2+dy**2
410 IF (ar.GE.(-armn) .AND. dsqi.GE.dsqmn) GO TO 230
411 jpmn = jp2
412 dxmn = dx
413 dymn = dy
414 dsqmn = dsqi
415 armn = dsqmn*ratio
416 end IF
417230 ar = dy*dxmx-dx*dymx
418 IF (ar .LT. (-armx)) cycle
419 dsqi = dx**2+dy**2
420 IF (ar.LE.armx .AND. dsqi.GE.dsqmx) cycle
421 jpmx = jp2
422 dxmx = dx
423 dymx = dy
424 dsqmx = dsqi
425 armx = dsqmx*ratio
426 end DO
427 IF (jpmx .LT. jpmn) jpmx = jpmx+nl0
428 nsh = jpmn-1
429 IF (nsh > 0) then
430 !
431 ! - SHIFTS (ROTATES) THE IPL ARRAY TO HAVE THE INVISIBLE BORDER
432 ! - LINE SEGMENTS CONTAINED IN THE FIRST PART OF THE IPL ARRAY.
433 !
434 nsht3 = nsh*3
435 DO jp2t3=3,nsht3,3
436 jp3t3 = jp2t3+nlt3
437 ipl(jp3t3-2) = ipl(jp2t3-2)
438 ipl(jp3t3-1) = ipl(jp2t3-1)
439 ipl(jp3t3) = ipl(jp2t3)
440 end DO
441 DO jp2t3=3,nlt3,3
442 jp3t3 = jp2t3+nsht3
443 ipl(jp2t3-2) = ipl(jp3t3-2)
444 ipl(jp2t3-1) = ipl(jp3t3-1)
445 ipl(jp2t3) = ipl(jp3t3)
446 end DO
447 jpmx = jpmx-nsh
448 !
449 ! - ADDS TRIANGLES TO THE IPT ARRAY, UPDATES BORDER LINE
450 ! - SEGMENTS IN THE IPL ARRAY, AND SETS FLAGS FOR THE BORDER
451 ! - LINE SEGMENTS TO BE REEXAMINED IN THE IWL ARRAY.
452 !
453 end IF
454
455 jwl = 0
456 l310 : DO jp2=jpmx,nl0
457 jp2t3 = jp2*3
458 ipl1 = ipl(jp2t3-2)
459 ipl2 = ipl(jp2t3-1)
460 it = ipl(jp2t3)
461 !
462 ! - - ADDS A TRIANGLE TO THE IPT ARRAY.
463 !
464 nt0 = nt0+1
465 ntt3 = ntt3+3
466 ipt(ntt3-2) = ipl2
467 ipt(ntt3-1) = ipl1
468 ipt(ntt3) = ip1
469 !
470 ! - - UPDATES BORDER LINE SEGMENTS IN THE IPL ARRAY.
471 !
472 IF (jp2 == jpmx) then
473 ipl(jp2t3-1) = ip1
474 ipl(jp2t3) = nt0
475 end IF
476 IF (jp2 == nl0) then
477 nln = jpmx+1
478 nlnt3 = nln*3
479 ipl(nlnt3-2) = ip1
480 ipl(nlnt3-1) = ipl(1)
481 ipl(nlnt3) = nt0
482 end IF
483 !
484 ! - - DETERMINES THE VERTEX THAT DOES NOT LIE ON THE BORDER
485 ! - - LINE SEGMENTS.
486 !
487 itt3 = it*3
488 ipti = ipt(itt3-2)
489 IF (ipti.NE.ipl1 .AND. ipti.NE.ipl2) GO TO 300
490 ipti = ipt(itt3-1)
491 IF (ipti.NE.ipl1 .AND. ipti.NE.ipl2) GO TO 300
492 ipti = ipt(itt3)
493 !
494 ! - - CHECKS IF THE EXCHANGE IS NECESSARY.
495 !
496300 IF (conxch_simc(co%X,co%Y,ip1,ipti,ipl1,ipl2) .EQ. 0) cycle l310
497 !
498 ! - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
499 !
500 ipt(itt3-2) = ipti
501 ipt(itt3-1) = ipl1
502 ipt(itt3) = ip1
503 ipt(ntt3-1) = ipti
504 IF (jp2 .EQ. jpmx) ipl(jp2t3) = it
505 IF (jp2.EQ.nl0 .AND. ipl(3).EQ.it) ipl(3) = nt0
506 !
507 ! - - SETS FLAGS IN THE IWL ARRAY.
508 !
509 jwl = jwl+4
510 iwl(jwl-3) = ipl1
511 iwl(jwl-2) = ipti
512 iwl(jwl-1) = ipti
513 iwl(jwl) = ipl2
514 end DO l310
515 nl0 = nln
516 nlt3 = nlnt3
517 nlf = jwl/2
518 IF (nlf .EQ. 0) cycle l400
519 !
520 ! - IMPROVES TRIANGULATION.
521 !
522 ntt3p3 = ntt3+3
523 DO irep=1,nrep
524 l370 : DO ilf=1,nlf
525 ilft2 = ilf*2
526 ipl1 = iwl(ilft2-1)
527 ipl2 = iwl(ilft2)
528 !
529 ! - - LOCATES IN THE IPT ARRAY TWO TRIANGLES ON BOTH SIDES OF
530 ! - - THE FLAGGED LINE SEGMENT.
531 !
532 ntf = 0
533 DO itt3r=3,ntt3,3
534 itt3 = ntt3p3-itt3r
535 ipt1 = ipt(itt3-2)
536 ipt2 = ipt(itt3-1)
537 ipt3 = ipt(itt3)
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
540 ntf = ntf+1
541 itf(ntf) = itt3/3
542 IF (ntf .EQ. 2) GO TO 330
543 end DO
544 IF (ntf .LT. 2) cycle l370
545 !
546 ! - - DETERMINES THE VERTEXES OF THE TRIANGLES THAT DO NOT LIE
547 ! - - ON THE LINE SEGMENT.
548 !
549330 it1t3 = itf(1)*3
550 ipti1 = ipt(it1t3-2)
551 IF (ipti1.NE.ipl1 .AND. ipti1.NE.ipl2) GO TO 340
552 ipti1 = ipt(it1t3-1)
553 IF (ipti1.NE.ipl1 .AND. ipti1.NE.ipl2) GO TO 340
554 ipti1 = ipt(it1t3)
555340 it2t3 = itf(2)*3
556 ipti2 = ipt(it2t3-2)
557 IF (ipti2.NE.ipl1 .AND. ipti2.NE.ipl2) GO TO 350
558 ipti2 = ipt(it2t3-1)
559 IF (ipti2.NE.ipl1 .AND. ipti2.NE.ipl2) GO TO 350
560 ipti2 = ipt(it2t3)
561 !
562 ! - - CHECKS IF THE EXCHANGE IS NECESSARY.
563 !
564350 IF (conxch_simc(co%X,co%Y,ipti1,ipti2,ipl1,ipl2) .EQ. 0) cycle l370
565 !
566 ! - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
567 !
568 ipt(it1t3-2) = ipti1
569 ipt(it1t3-1) = ipti2
570 ipt(it1t3) = ipl1
571 ipt(it2t3-2) = ipti2
572 ipt(it2t3-1) = ipti1
573 ipt(it2t3) = ipl2
574 !
575 ! - - SETS NEW FLAGS.
576 !
577 jwl = jwl+8
578 iwl(jwl-7) = ipl1
579 iwl(jwl-6) = ipti1
580 iwl(jwl-5) = ipti1
581 iwl(jwl-4) = ipl2
582 iwl(jwl-3) = ipl2
583 iwl(jwl-2) = ipti2
584 iwl(jwl-1) = ipti2
585 iwl(jwl) = ipl1
586 DO jlt3=3,nlt3,3
587 iplj1 = ipl(jlt3-2)
588 iplj2 = ipl(jlt3-1)
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)
593 end DO
594 end DO l370
595 nlfc = nlf
596 nlf = jwl/2
597 IF (nlf .EQ. nlfc) cycle l400
598 !
599 ! - - RESETS THE IWL ARRAY FOR THE NEXT ROUND.
600 !
601 jwl = 0
602 jwl1mn = (nlfc+1)*2
603 nlft2 = nlf*2
604 DO jwl1=jwl1mn,nlft2,2
605 jwl = jwl+2
606 iwl(jwl-1) = iwl(jwl1-1)
607 iwl(jwl) = iwl(jwl1)
608 end DO
609 nlf = jwl/2
610 end DO
611end DO l400
612
613call l4f_log(l4f_debug,"end triangulation appending")
614
615 !
616 ! REARRANGE THE IPT ARRAY SO THAT THE VERTEXES OF EACH TRIANGLE
617 ! ARE LISTED COUNTER-CLOCKWISE.
618 !
619DO itt3=3,ntt3,3
620 ip1 = ipt(itt3-2)
621 ip2 = ipt(itt3-1)
622 ip3 = ipt(itt3)
623 IF (side(co(ip1),co(ip2),co(ip3)) .GE. 10.0) cycle
624 ipt(itt3-2) = ip2
625 ipt(itt3-1) = ip1
626end DO
627call l4f_log(l4f_debug,"end triangulation rearranging")
628
629nt = nt0
630nl = nl0
631
632call l4f_log(l4f_debug,"end triangulation")
633
634RETURN
635
636contains
637
638!!$double precision function DSQF(U1,V1,U2,V2)
639!!$double precision,intent(in) :: U1,V1,U2,V2
640!!$
641!!$DSQF = (U2-U1)**2+(V2-V1)**2
642!!$end function DSQF
643
644
645elemental double precision function vdsqf(co1,co2)
646type(xy),intent(in) :: co1,co2
647
648vdsqf = (co2%x-co1%x)**2+(co2%y-co1%y)**2
649if (vdsqf == 0.d0) vdsqf = huge(vdsqf)
650end function vdsqf
651
652
653!!$double precision function SIDE(U1,V1,U2,V2,U3,V3)
654!!$double precision,intent(in):: U1,V1,U2,V2,U3,V3
655!!$
656!!$SIDE = (V3-V1)*(U2-U1)-(U3-U1)*(V2-V1)
657!!$end function SIDE
658
659double precision function side(co1,co2,co3)
660type(xy),intent(in):: co1,co2,co3
661
662side = (co3%y-co1%y)*(co2%x-co1%x)-(co3%x-co1%x)*(co2%y-co1%y)
663end function side
664
665end function contng_simc
666
667
673INTEGER FUNCTION conxch_simc (X,Y,I1,I2,I3,I4)
674
675double precision,intent(in) :: x(:), y(:)
678integer,intent(in) :: i1,i2,i3,i4
679
680double precision :: x0(4), y0(4)
681double precision :: c2sq,c1sq,a3sq,b2sq,b3sq,a1sq,a4sq,b1sq,b4sq,a2sq,c4sq,c3sq
682integer :: idx
683double precision :: s1sq,s2sq,s3sq,s4sq,u1,u2,u3,u4
684
685equivalence(c2sq,c1sq),(a3sq,b2sq),(b3sq,a1sq),(a4sq,b1sq), &
686 (b4sq,a2sq),(c4sq,c3sq)
687
688x0(1) = x(i1)
689y0(1) = y(i1)
690x0(2) = x(i2)
691y0(2) = y(i2)
692x0(3) = x(i3)
693y0(3) = y(i3)
694x0(4) = x(i4)
695y0(4) = y(i4)
696idx = 0
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))
699IF (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
713end IF
714conxch_simc = idx
715RETURN
716
717END FUNCTION conxch_simc
718
719end 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.