libsim  Versione 7.2.6
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 
60 module space_utilities
61 use log4fortran
64 implicit none
65 
66 type :: triangles
67  integer,pointer :: ipt(:) => null(), ipl(:) => null()
68  integer :: nt=imiss,nl=imiss
69 end type triangles
70 
71 type :: xy
72  double precision :: x,y
73 end type xy
74 
77 INTERFACE delete
78  MODULE PROCEDURE triangles_delete
79 END INTERFACE
80 
81 INTERFACE triangles_compute
82  MODULE PROCEDURE triangles_compute_r, triangles_compute_d, triangles_compute_c
83 END INTERFACE
84 
85 private
86 public triangles, triangles_new, delete, triangles_compute, xy
87 
88 
89 contains
90 
92 function triangles_new(ndp) result(this)
93 type(triangles) :: this
94 integer,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 
101 if (c_e(ndp) .and. ndp >= 3)then
102  allocate(this%ipt(6*ndp-15), this%ipl(6*ndp))
103 else
104  this%nt=0
105  this%nl=0
106 end if
107 return
108 end function triangles_new
109 
110 
112 subroutine triangles_delete(this)
113 type(triangles) :: this
114 
115 if (associated(this%ipt)) deallocate(this%ipt)
116 if (associated(this%ipl)) deallocate(this%ipl)
117 
118 this%nt=imiss
119 this%nl=imiss
120 
121 end subroutine triangles_delete
122 
123 
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))
129 
130 if (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)
134 end if
135 end function triangles_compute_r
136 
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))
142 
143 if (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)
147 end if
148 end function triangles_compute_d
149 
150 integer function triangles_compute_c (co,tri)
151 type (xy),intent(in) :: co(:)
152 type (triangles),intent(inout) :: tri
153 
154 if (tri%nt /= 0) then
155  triangles_compute_c = contng_simc(co,tri%NT,tri%IPT,tri%NL,tri%IPL)
156 end if
157 end function triangles_compute_c
158 
159 
173 integer function contng_simc (co,NT,IPT,NL,IPL)
174 
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(:)
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.
197 integer :: IWL(18*size(co)), IWP(size(co))
198 
199 double precision :: WK(size(co))
200 integer :: ITF(2), NREP=100
201 real :: RATIO=1.0e-6
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
209 logical :: err
210 
211 integer ::i,mloc(size(co)-1),mlocall(1),mlocv(1)
212 type(xy) :: dmp
213 
214  !
215  ! PRELIMINARY PROCESSING
216  !
217 
218 nt = imiss
219 nl = imiss
220 
221 contng_simc=0
222 ndp=size(co)
223 ndpm1 = ndp-1
224  !
225  ! DETERMINES THE CLOSEST PAIR OF DATA POINTS AND THEIR MIDPOINT.
226  !
227 
228 call l4f_log(l4f_debug,"start triangulation")
229 
230 do i=1,size(co)-1
231  mlocv=minloc(vdsqf(co(i),co(i+1:)))+i
232  mloc(i)=mlocv(1)
233 end do
234 
235 mlocall=minloc((/(vdsqf(co(i),co(mloc(i))),i=1,size(mloc))/))
236 
237 dsqmn = vdsqf(co(mlocall(1)),co(mloc(mlocall(1))))
238 ipmn1 = mlocall(1)
239 ipmn2 = mloc(mlocall(1))
240 
241 call l4f_log(l4f_debug,"end triangulation closest pair")
242 !!$print *, DSQMN, IPMN1, IPMN2
243 
244 IF (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
251 end 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 
284 dsq12 = dsqmn
285 dmp%x = (co(ipmn1)%x+co(ipmn2)%x)/2.0
286 dmp%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  !
292 jp1 = 2
293 DO 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))
298 end DO
299 DO 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)
311 end DO
312 
313 call 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  !
319 ar = dsq12*ratio
320 x1 = co(ipmn1)%x
321 y1 = co(ipmn1)%y
322 dx21 = co(ipmn2)%x-x1
323 dy21 = co(ipmn2)%y-y1
324 
325 err=.true.
326 DO 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
332 end DO
333 if (err) then
334  call l4f_log(l4f_debug,"CONTNG - ALL COLLINEAR DATA POINTS")
335  contng_simc=2
336  return
337 end if
338 IF (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
346 end IF
347 call 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  !
355 ip1 = ipmn1
356 ip2 = ipmn2
357 ip3 = iwp(3)
358 IF (side(co(ip1),co(ip2),co(ip3)) < 10.0) then
359  ip1 = ipmn2
360  ip2 = ipmn1
361 end IF
362 
363 nt0 = 1
364 ntt3 = 3
365 ipt(1) = ip1
366 ipt(2) = ip2
367 ipt(3) = ip3
368 nl0 = 3
369 nlt3 = 9
370 ipl(1) = ip1
371 ipl(2) = ip2
372 ipl(3) = 1
373 ipl(4) = ip2
374 ipl(5) = ip3
375 ipl(6) = 1
376 ipl(7) = ip3
377 ipl(8) = ip1
378 ipl(9) = 1
379 
380 call l4f_log(l4f_debug,"end triangulation first triangle")
381 
382  !
383  ! ADDS THE REMAINING (NDP-3) DATA POINTS, ONE BY ONE.
384  !
385 l400 : 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
417 230 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  !
496 300 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  !
549 330 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)
555 340 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  !
564 350 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
611 end DO l400
612 
613 call 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  !
619 DO 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
626 end DO
627 call l4f_log(l4f_debug,"end triangulation rearranging")
628 
629 nt = nt0
630 nl = nl0
631 
632 call l4f_log(l4f_debug,"end triangulation")
633 
634 RETURN
635 
636 contains
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 
645 elemental double precision function vdsqf(co1,co2)
646 type(xy),intent(in) :: co1,co2
647 
648 vdsqf = (co2%x-co1%x)**2+(co2%y-co1%y)**2
649 if (vdsqf == 0.d0) vdsqf = huge(vdsqf)
650 end 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 
659 double precision function side(co1,co2,co3)
660 type(xy),intent(in):: co1,co2,co3
661 
662 side = (co3%y-co1%y)*(co2%x-co1%x)-(co3%x-co1%x)*(co2%y-co1%y)
663 end function side
664 
665 end function contng_simc
666 
667 
673 INTEGER FUNCTION conxch_simc (X,Y,I1,I2,I3,I4)
674 
675 double precision,intent(in) :: x(:), y(:)
678 integer,intent(in) :: i1,i2,i3,i4
679 
680 double precision :: x0(4), y0(4)
681 double precision :: c2sq,c1sq,a3sq,b2sq,b3sq,a1sq,a4sq,b1sq,b4sq,a2sq,c4sq,c3sq
682 integer :: idx
683 double precision :: s1sq,s2sq,s3sq,s4sq,u1,u2,u3,u4
684 
685 equivalence(c2sq,c1sq),(a3sq,b2sq),(b3sq,a1sq),(a4sq,b1sq), &
686  (b4sq,a2sq),(c4sq,c3sq)
687 
688 x0(1) = x(i1)
689 y0(1) = y(i1)
690 x0(2) = x(i2)
691 y0(2) = y(i2)
692 x0(3) = x(i3)
693 y0(3) = y(i3)
694 x0(4) = x(i4)
695 y0(4) = y(i4)
696 idx = 0
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
713 end IF
714 conxch_simc = idx
715 RETURN
716 
717 END FUNCTION conxch_simc
718 
719 end 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.