!npox,ntrix

! ================================================
!   written by MATTIA DE MICHIELI VITTURI 2002
! ================================================
!       last modified 19/11/2003
! ================================================
SUBROUTINE triangola(npox,ntrix,xRR,yRR,nodiest,noditot,triang,nt)
  USE fund_const
  IMPLICIT none
  INTEGER, INTENT(IN) :: npox,ntrix ! max number of points and triangles
  INTEGER, INTENT(IN) :: nodiest ! number of nodes on the boundary
  INTEGER, INTENT(OUT) :: noditot ! number of nodes
  DOUBLE PRECISION, DIMENSION(npox), INTENT(INOUT) :: xRR,yRR 
! input: nodiest boundary nodes
! output: all noditot nodes
  INTEGER,INTENT(OUT) :: nt ! number of triangles
  INTEGER, INTENT(OUT) :: triang(ntrix,4)
! END INTERFACE
! =========================================================================
  DOUBLE PRECISION xR(0:npox),yR(0:npox)
  INTEGER circmax,seg(ntrix,3,2),triangnew(4),nodo(4)
  INTEGER segnew3(3,2),segmento(npox,npox),nodoopposto(3)
  INTEGER escimax,h,l,i,m,j
  INTEGER intersez,badnodi,cambia(ntrix),vicini(3),quantivicini
  INTEGER quanti2,triangmax(4)
  DOUBLE PRECISION raggio,raggiont,xmax1,xmax2,xmax3,ymax1,ymax2,ymax3
  DOUBLE PRECISION angolo,raggiocirccentr(ntrix)
  DOUBLE PRECISION a11,a12,a21,a22
  DOUBLE PRECISION xcirccentr(ntrix),yverif1,yverif2,yverif3
  DOUBLE PRECISION ycirccentr(ntrix),xverif1,xverif2,xverif3
  DOUBLE PRECISION xbaric(ntrix),ybaric(ntrix)
  DOUBLE PRECISION density(npox),distmin
  DOUBLE PRECISION ITE, TA(2)
  DOUBLE PRECISION xcentro,ycentro,raggiomax,area2,etime
  INTEGER iun
  DOUBLE PRECISION cos1,cos2,sen1,sen2,angolo1,angolo2
  INTEGER exterior
  INTEGER :: count1,count2 ! counters
  DOUBLE PRECISION, PARAMETER :: raggiofix=0.58d0 
  DOUBLE PRECISION, PARAMETER :: epsilon0=1.D-3
  DOUBLE PRECISION, PARAMETER :: eps1=1.D-3
  INTEGER, PARAMETER :: ntmax=5000
  
  DO i = 1,nodiest
     xR(i)=xRR(i)
     yR(i)=yRR(i)
  ENDDO
  xR(nodiest+1)=xR(1)
  yR(nodiest+1)=yR(1)
  xR(0)=xR(nodiest)
  yR(0)=yR(nodiest)
  
  noditot=nodiest
      
  raggiocirccentr=0.d0      
            
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&      
! write on file external nodes
  CALL filopn(iun,'bordo_esterno','unknown')
  DO j =1, nodiest
!         WRITE(iun,199)j,xR(j),yR(j)
     WRITE(iun,*)j,xR(j),yR(j)
  ENDDO

199 FORMAT(i5,1x,f20.17,1x,f20.17)
  CALL filclo(iun,' ')
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&      


!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&      
! INITIAL TRIANGULATION
  nt=0
  DO i=1,noditot
     DO j=1,noditot
        segmento(i,j)=0
        IF ((j-i.EQ.1).OR.(j-i.EQ.nodiest-1)) segmento(i,j)=1
     ENDDO
  ENDDO
  
10 DO i=1,noditot-1
     DO j=i+1,noditot
        IF (segmento(i,j).EQ.1) THEN 
           GOTO 5
        ENDIF
     ENDDO
  ENDDO
  
  GOTO 39
  
! THEN we must add a triangle using the two nodes i,j     
5 nt=nt+1
  IF(nt.GT.ntrix)THEN
     write(*,*)'triangola: nt too large ',nt,ntrix
     STOP
  ENDIF
  
  escimax=noditot+1
  DO h=j+1,noditot
     exterior=0
     IF ((h.ne.j).AND.(h.ne.i)) THEN 
        
        nodo(1)=i
        nodo(2)=j
        nodo(3)=h
        nodo(4)=i
        
        DO m=1,3
           a11=yR(nodo(m)-1)-yR(nodo(m))
           a12=yR(nodo(m)+1)-yR(nodo(m))
           a21=xR(nodo(m)-1)-xR(nodo(m))
           a22=xR(nodo(m)+1)-xR(nodo(m))
           sen1=(a11*a22-a12*a21)/(sqrt(a11**2+a21**2)*sqrt(a12**2+a22**2))
           cos1=(a11*a12+a21*a22)/(sqrt(a11**2+a21**2)*sqrt(a12**2+a22**2))
           angolo1=mod(atan2(sen1,cos1)+2.d0*pig,2.d0*pig)
           a11=yR(nodo(m)-1)-yR(nodo(m))
           a12=yR(nodo(m+1))-yR(nodo(m))
           a21=xR(nodo(m)-1)-xR(nodo(m))
           a22=xR(nodo(m+1))-xR(nodo(m))
           sen2=(a11*a22-a12*a21)/(sqrt(a11**2+a21**2)*sqrt(a12**2+a22**2))
           cos2=(a11*a12+a21*a22)/(sqrt(a11**2+a21**2)*sqrt(a12**2+a22**2))
           angolo2=mod(atan2(sen2,cos2)+2.d0*pig,2.d0*pig)
           IF (angolo2.GT.angolo1+1.d-7) THEN
              exterior=1
           ENDIF
        ENDDO
        a11=yR(j)-yR(i)
        a12=-(yR(h)-yR(i))
        a21=-(xR(j)-xR(i))
        a22=xR(h)-xR(i)
        IF (abs(a11*a22-a12*a21)/(sqrt(a11**2+a21**2)&
             &  *sqrt(a12**2+a22**2)).le.epsilon0) exterior=2
        
        IF (exterior.EQ.0) THEN
           triangnew(1)=i
           triangnew(2)=j
           triangnew(3)=h
           triangnew(4)=i
           
           CALL nuovotriangolo2(triangnew,segnew3,&
                & xR(i),yR(i),xR(j),yR(j),xR(h),yR(h),&
                & xcentro,ycentro,raggiont)
           raggio=raggiont
           
           CALL verificaint(npox,ntrix,triang,nt,xR,yR,intersez,xcentro,&
                & ycentro,raggio,badnodi,noditot,triangnew,segmento)
           
           IF ((intersez.EQ.0).AND.((badnodi.LT.escimax)&
                & .OR.((badnodi.EQ.escimax).AND. &
                & (raggiont.lt.raggiocirccentr(nt)*(1.d0+ &
                & epsilon(1.d0)*100.d0))))) THEN
              escimax=badnodi
              
              CALL nuovotriangolo(triangnew,segnew3,&
                   & xR(i),yR(i),xR(j),yR(j),xR(h),yR(h),&
                   & xcentro,ycentro,raggiont,&
                   & xbaric(nt),ybaric(nt))
              raggio=raggiont
              
              triang(nt,1)=triangnew(1)
              triang(nt,2)=triangnew(2)
              triang(nt,3)=triangnew(3)
              xcirccentr(nt)=xcentro
              ycirccentr(nt)=ycentro
              raggiocirccentr(nt)=raggiont
              DO l=1,3
                 seg(nt,l,1)=segnew3(l,1)
                 seg(nt,l,2)=segnew3(l,2)
              ENDDO
           ENDIF
        ENDIF
     ENDIF
     
  ENDDO
!      write(*,*)'nt=',nt,triang(nt,1),triang(nt,2),triang(nt,3)   
  IF(triang(nt,1).EQ.0) THEN
     WRITE(*,*)'triangola: error for segment',i,j,&
          &' triang(nt,1)=',triang(nt,1)
     STOP
  ENDIF
  segmento(triang(nt,1),triang(nt,2))=segmento(triang(nt,1),triang(nt,2))+1
  segmento(triang(nt,2),triang(nt,3))=segmento(triang(nt,2),triang(nt,3))+1
  segmento(triang(nt,1),triang(nt,3))=segmento(triang(nt,1),triang(nt,3))+1
  
  GOTO 10
  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&      
  
  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ! OPTIMIZATION OF THE INITIAL TRIANGULATION 
  ! (NEEDED BECAUSE DOMAIN NOT CONVEX)
39 count2 = 0
  
40 count2 = count2+1
  if(count2.GT.2**nt) THEN
     GOTO 36
  ENDIF
  
  DO i=1,nt
     
     CALL cercatriangolivicini(npox,ntrix,triang,nt,i,quantivicini,  &
          &   vicini,seg,nodoopposto)
     
     DO j=1,quantivicini
        
        IF (sqrt((xR(nodoopposto(j))-xcirccentr(i))           &
             &  **2+(yR(nodoopposto(j))-ycirccentr(i))**2)    &
             &  .LT.raggiocirccentr(i)*(1.d0-eps1)) THEN
           
           CALL listaswap(npox,ntrix,triang,xR,yR,xcirccentr,            &
                &   seg,ycirccentr,raggiocirccentr,           &
                &   i,quanti2,cambia,nt,eps1)
           
           CALL swappalista(npox,ntrix,triang,xR,yR,xcirccentr,seg,      &
                &   ycirccentr,raggiocirccentr,quanti2,cambia,&
                &   xbaric,ybaric,eps1)
           GOTO 40
           
        ENDIF
        
     ENDDO
     
  ENDDO

  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ! COMPUTE DENSITY ON THE BOUNDARY
36 DO i=1,noditot
     
     density(i)=0.5d0*sqrt((xR(i)-xR(i-1))**2+(yR(i)-yR(i-1))**2)+ &
          &   0.5d0*sqrt((xR(i)-xR(i+1))**2+(yR(i)-yR(i+1))**2)
     
  ENDDO  
  
  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ! STARTING TRIANGULATION REFINEMENT
  raggiomax=2.d0*raggiofix
  
  !      DO while (raggiomax.ge.raggiofix)
  DO while ((raggiomax.ge.raggiofix).AND.(nt.LT.ntmax-2))
     
     raggiomax=0
     DO i=1,nt
        xverif1=xR(triang(i,1))
        yverif1=yR(triang(i,1))
        xverif2=xR(triang(i,2))
        yverif2=yR(triang(i,2))
        xverif3=xR(triang(i,3))
        yverif3=yR(triang(i,3))
        
        
        distmin=1.d5
        DO j=1,3
           distmin = min(distmin,sqrt((xbaric(i)-xR(triang(i,j)))**2 &
                &  +(ybaric(i)-yR(triang(i,j)))**2))
        ENDDO
        
        area2 = 3.d0*distmin/( density(triang(i,1)) + &
             & density(triang(i,2))+ density(triang(i,3)) )
        
        IF (area2.ge.raggiomax) THEN
           raggiomax=area2
           circmax=i
        ENDIF
     ENDDO
     
     noditot=noditot+1
     
     xR(noditot) = xbaric(circmax)
     yR(noditot) = ybaric(circmax)
     
     density(noditot)=(density(triang(circmax,1))              &
          &  +density(triang(circmax,2))+density(triang(circmax,3)))/3
     
     xmax1=xR(triang(circmax,1))
     ymax1=yR(triang(circmax,1))
     xmax2=xR(triang(circmax,2))
     ymax2=yR(triang(circmax,2))
     xmax3=xR(triang(circmax,3))
     ymax3=yR(triang(circmax,3))
     triangmax(1)=triang(circmax,1)
     triangmax(2)=triang(circmax,2)
     triangmax(3)=triang(circmax,3)
     triangmax(4)=triang(circmax,1)
     
     DO i=circmax,nt-1
        triang(i,1)=triang(i+1,1)
        triang(i,2)=triang(i+1,2)
        triang(i,3)=triang(i+1,3)
        triang(i,4)=triang(i+1,3)
        seg(i,1,1)=seg(i+1,1,1)
        seg(i,1,2)=seg(i+1,1,2)
        seg(i,2,1)=seg(i+1,2,1)
        seg(i,2,2)=seg(i+1,2,2)
        seg(i,3,1)=seg(i+1,3,1)
        seg(i,3,2)=seg(i+1,3,2)
        raggiocirccentr(i)=raggiocirccentr(i+1)
        xcirccentr(i)=xcirccentr(i+1)
        ycirccentr(i)=ycirccentr(i+1)
        xbaric(i)=xbaric(i+1)
        ybaric(i)=ybaric(i+1)
        
     ENDDO
     nt=nt-1
     DO i=1,3
        nt=nt+1
        triangnew(1)=triangmax(i)
        triangnew(2)=triangmax(i+1)
        triangnew(3)=noditot
        triangnew(4)=triangmax(i)
        CALL nuovotriangolo(triangnew,segnew3,      &
             &         xR(triangnew(1)),yR(triangnew(1)),       &
             &         xR(triangnew(2)),yR(triangnew(2)),       &
             &         xR(triangnew(3)),yR(triangnew(3)),       &
             &         xcentro,ycentro,raggiont,                &
             &         xbaric(nt),ybaric(nt))
        
        DO l=1,3
           seg(nt,l,1)=segnew3(l,1)
           seg(nt,l,2)=segnew3(l,2)
        ENDDO
        triang(nt,1)=triangnew(1)
        triang(nt,2)=triangnew(2)
        triang(nt,3)=triangnew(3)
        xcirccentr(nt)=xcentro
        ycirccentr(nt)=ycentro
        raggiocirccentr(nt)=raggiont
        CALL cercatriangolivicini(npox,ntrix,triang,nt,nt,quantivicini, &
             &  vicini,seg,nodoopposto)
        DO j=1,quantivicini
           IF (sqrt((xR(nodoopposto(j))-xcirccentr(nt))**2  &
                &  +(yR(nodoopposto(j))-ycirccentr(nt))**2)   &
                &  .LT.raggiocirccentr(nt)*(1.d0-eps1)) THEN
              CALL listaswap(npox,ntrix,triang,xR,yR,xcirccentr, &
                   &   seg,ycirccentr,raggiocirccentr, &
                   &   nt,quanti2,cambia,nt,eps1)
              CALL swappalista(npox,ntrix,triang,xR,yR,xcirccentr,seg, &
                   &   ycirccentr,raggiocirccentr,quanti2,cambia, &
                   &   xbaric,ybaric,eps1)
           ENDIF
        ENDDO
     ENDDO
  ENDDO
  CALL smoothing(npox,ntrix,eps1,nodiest,noditot,nt,seg,xR,yR,density,triang)
  DO i = 1,noditot
     xRR(i)=xR(i)
     yRR(i)=yR(i)
  ENDDO
  
  RETURN
END SUBROUTINE triangola
! &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!                       !!!!!!!!!!!
!!!!!!!!!          NUOVOTRIANGOLO   !!!!!!!
!!!!!!!!!                           !!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE nuovotriangolo(triangnew1,segnew,xR1,yR1,xR2,yR2, &
     &          xR3,yR3,xcentronew,ycentronew,raggionew,xbaricx,ybaricx)
  IMPLICIT NONE
  INTEGER, INTENT(INOUT) :: triangnew1(4)
  INTEGER, INTENT(OUT) :: segnew(3,2)
  DOUBLE PRECISION, INTENT(IN) :: xR1,yR1,xR2,yR2,xR3,yR3
  DOUBLE PRECISION, INTENT(OUT)  :: xcentronew,ycentronew,raggionew
  DOUBLE PRECISION, INTENT(OUT)  :: xbaricx,ybaricx
  ! END INTERFACE
  DOUBLE PRECISION :: t,s
  DOUBLE PRECISION :: a11,a12,a21,a22,b1,b2
  DOUBLE PRECISION :: xcirccentr1,xcirccentr2,ycirccentr1,ycirccentr2
  ! =========================================================================
  
  IF (triangnew1(1).GT.triangnew1(2)) THEN
     triangnew1(4)=triangnew1(1)
     triangnew1(1)=triangnew1(2)
     triangnew1(2)=triangnew1(4)
  ENDIF
  IF (triangnew1(1).GT.triangnew1(3)) THEN
     triangnew1(4)=triangnew1(1)
     triangnew1(1)=triangnew1(3)
     triangnew1(3)=triangnew1(4)
  ENDIF
  IF (triangnew1(2).GT.triangnew1(3)) THEN
     triangnew1(4)=triangnew1(2)
     triangnew1(2)=triangnew1(3)
     triangnew1(3)=triangnew1(4)
  ENDIF
  triangnew1(4)=triangnew1(1)
  segnew(1,1)=triangnew1(1)
  segnew(1,2)=triangnew1(2)
  segnew(2,1)=triangnew1(2)
  segnew(2,2)=triangnew1(3) 
  segnew(3,1)=triangnew1(1)
  segnew(3,2)=triangnew1(3)
  a11=yR2-yR1
  a12=-(yR3-yR1)
  a21=-(xR2-xR1)
  a22=xR3-xR1
  b1=0.5*(xR3-xR2)
  b2=0.5*(yR3-yR2)
  
  IF (a11*a22-a12*a21.EQ.0.d0) THEN
     write(*,*) 'nuovotriangolo: zero denominator'
     stop
  ENDIF
  
  t=(b1*a22-a12*b2)/(a11*a22-a12*a21)
  s=(a11*b2-b1*a21)/(a11*a22-a12*a21)
  xcirccentr1=(xR1+xR2)/2.d0+t*(yR2-yR1)
  ycirccentr1=(yR1+yR2)/2.d0-t*(xR2-xR1)
  xcirccentr2=(xR1+xR3)/2.d0+s*(yR3-yR1)
  ycirccentr2=(yR1+yR3)/2.d0-s*(xR3-xR1)
  xcentronew=0.5d0*(xcirccentr1+xcirccentr2)
  ycentronew=0.5d0*(ycirccentr1+ycirccentr2)
  raggionew=sqrt((xcentronew-xR1)**2+(ycentronew-yR1)**2)
  xbaricx=(xR1+xR2+xR3)/3.d0
  ybaricx=(yR1+yR2+yR3)/3.d0
  
  RETURN
END SUBROUTINE nuovotriangolo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!                       !!!!!!!!!!!
!!!!!!!!!     NUOVOTRIANGOLO 2       !!!!!!!!!!!
!!!!!!!!!                         !!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE nuovotriangolo2(triangnew1,segnew,xR1,yR1,xR2,yR2, &
     &    xR3,yR3,xcentronew,ycentronew,raggionew)
  IMPLICIT NONE
  INTEGER, INTENT(INOUT) :: triangnew1(4)
  INTEGER, INTENT(OUT) :: segnew(3,2)
  DOUBLE PRECISION, INTENT(IN) :: xR1,yR1,xR2,yR2,xR3,yR3
  DOUBLE PRECISION, INTENT(OUT) :: xcentronew,ycentronew,raggionew
! END INTERFACE
  DOUBLE PRECISION :: t,s
  DOUBLE PRECISION :: a11,a12,a21,a22,b1,b2
  DOUBLE PRECISION :: xcirccentr1,xcirccentr2,ycirccentr1,ycirccentr2
  ! =========================================================================
  
  
  IF (triangnew1(1).GT.triangnew1(2)) THEN
     triangnew1(4)=triangnew1(1)
     triangnew1(1)=triangnew1(2)
     triangnew1(2)=triangnew1(4)
  ENDIF
  IF (triangnew1(1).GT.triangnew1(3)) THEN
     triangnew1(4)=triangnew1(1)
     triangnew1(1)=triangnew1(3)
     triangnew1(3)=triangnew1(4)
  ENDIF
  IF (triangnew1(2).GT.triangnew1(3)) THEN
     triangnew1(4)=triangnew1(2)
     triangnew1(2)=triangnew1(3)
     triangnew1(3)=triangnew1(4)
  ENDIF
  triangnew1(4)=triangnew1(1)
  segnew(1,1)=triangnew1(1)
  segnew(1,2)=triangnew1(2)
  segnew(2,1)=triangnew1(2)
  segnew(2,2)=triangnew1(3) 
  segnew(3,1)=triangnew1(1)
  segnew(3,2)=triangnew1(3)
  a11=yR2-yR1
  a12=-(yR3-yR1)
  a21=-(xR2-xR1)
  a22=xR3-xR1
  b1=0.5*(xR3-xR2)
  b2=0.5*(yR3-yR2)
  
  IF (a11*a22-a12*a21.EQ.0) THEN
     write(*,*) 'nuovotriangolo2: zero denominator'
     stop
  ENDIF
  
  t=(b1*a22-a12*b2)/(a11*a22-a12*a21)
  s=(a11*b2-b1*a21)/(a11*a22-a12*a21)
  xcirccentr1=(xR1+xR2)/2+t*(yR2-yR1)
  ycirccentr1=(yR1+yR2)/2-t*(xR2-xR1)
  xcirccentr2=(xR1+xR3)/2+s*(yR3-yR1)
  ycirccentr2=(yR1+yR3)/2-s*(xR3-xR1)
  xcentronew=0.5*(xcirccentr1+xcirccentr2)
  ycentronew=0.5*(ycirccentr1+ycirccentr2)
  raggionew=sqrt((xcentronew-xR1)**2+(ycentronew-yR1)**2)
  RETURN
END SUBROUTINE nuovotriangolo2


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!                       !!!!!!!!!!!
!!!!!!!!!      LISTASWAP            !!!!!!!
!!!!!!!!!                       !!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE listaswap(npox,ntrix,triangnew,xRnew,yRnew,xcirccentrnew,seg, &
     &    ycirccentrnew,raggiocirccentrnew,primo,quanti2,cambia,nt,eps1)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: npox,ntrix ! max number of points and triangles
  INTEGER, INTENT(IN) :: triangnew(ntrix,4),primo
  INTEGER, INTENT(OUT):: quanti2,cambia(ntrix)
  INTEGER, INTENT(IN) :: nt,seg(ntrix,3,2)
  DOUBLE PRECISION, INTENT(IN) :: xRnew(0:npox),yRnew(0:npox)
  DOUBLE PRECISION, INTENT(IN) :: xcirccentrnew(ntrix)
  DOUBLE PRECISION, INTENT(IN) :: raggiocirccentrnew(ntrix),&
       & ycirccentrnew(ntrix)
  DOUBLE PRECISION, INTENT(IN) :: eps1
  ! END INTERFACE
  INTEGER :: i,j ! loop indexes
  INTEGER :: esiste,quantivicini,vicini(3),nodoopposto(3)
! =========================================================================    
  quanti2=1
  cambia(1)=primo
900 CALL cercatriangolivicini(npox,ntrix,triangnew,nt,cambia(1),quantivicini, &
         &    vicini,seg,nodoopposto)
  cambia(quanti2+1)=cambia(1)
  DO j=1,quanti2
     cambia(j)=cambia(j+1)
  ENDDO
  DO i=1,quantivicini
     esiste=0
     DO j=1,quanti2
        IF (vicini(i).EQ.cambia(j)) esiste=1
     ENDDO
     IF ((sqrt((xRnew(nodoopposto(i))-xcirccentrnew(primo))**2 &
          & +(yRnew(nodoopposto(i))-ycirccentrnew(primo))**2) &
          & .LT.raggiocirccentrnew(primo)*(1.d0-eps1)).AND.(esiste.EQ.0))THEN
        quanti2=quanti2+1
        DO j=2,quanti2
           cambia(quanti2+2-j)=cambia(quanti2+1-j)
        ENDDO
        cambia(1)=vicini(i)
     ENDIF
  ENDDO
  IF (cambia(1).ne.primo) GOTO 900
  RETURN
END SUBROUTINE listaswap


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!                       !!!!!!!!!!!
!!!!!!!!!      SWAPPALISTA           !!!!!!
!!!!!!!!!                       !!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     

SUBROUTINE swappalista(npox,ntrix,triangnew,xRnew,yRnew,xcirccentrnew,seg, &
     &    ycirccentrnew,raggiocirccentrnew,quanti,cambiain,xbaric,ybaric, &
     &    eps1)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: npox,ntrix ! max number of points and triangles
  INTEGER, INTENT(INOUT) :: triangnew(ntrix,4)
  INTEGER, INTENT(IN) :: quanti,cambiain(ntrix)
  INTEGER, INTENT(INOUT) :: seg(ntrix,3,2)
  DOUBLE PRECISION, INTENT(IN) :: xRnew(0:npox),yRnew(0:npox)
  DOUBLE PRECISION,INTENT(INOUT) :: xcirccentrnew(ntrix)
  DOUBLE PRECISION,INTENT(INOUT):: raggiocirccentrnew(ntrix),&
       & ycirccentrnew(ntrix)
  DOUBLE PRECISION,INTENT(INOUT):: xbaric(ntrix),ybaric(ntrix)
  DOUBLE PRECISION,INTENT(IN) :: eps1
  ! END INTERFACE
  INTEGER :: cambia(ntrix)
  INTEGER :: h,i,j,l ! loop indexes
  INTEGER :: count,ntg ! counters
  INTEGER :: nodo2,nodo4
  INTEGER :: segnew4(3,2),segnew5(3,2),primo,triangnew4(4)
  DOUBLE PRECISION :: xcirco,ycirco,raggiocirco
! =========================================================================
  
  count = 0
  DO h = 1,quanti
     cambia(h) = cambiain(h)   
  ENDDO
  
  ntg=0
800 primo=cambia(1)
  
  count = count+1
  if(count.GT.2**quanti) THEN
     !       write(*,*) 'there is a loop'
     GOTO 313
  ENDIF
  
  IF (ntg.LT.quanti) THEN
     DO i=2,quanti
        DO j=1,3
           DO h=1,3
              IF ((seg(primo,j,1).EQ.seg(cambia(i),h,1)).AND.&
                   &  (seg(primo,j,2).EQ.seg(cambia(i),h,2))) THEN
                 IF (h.EQ.1) THEN
                    nodo4=triangnew(cambia(i),3)
                 ELSEIF (h.EQ.2) THEN
                    nodo4=triangnew(cambia(i),1)
                 ELSE
                    nodo4=triangnew(cambia(i),2)
                 ENDIF
                 IF (j.EQ.1) THEN
                    nodo2=triangnew(primo,3)
                 ELSEIF (j.EQ.2) THEN
                    nodo2=triangnew(primo,1)
                 ELSE
                    nodo2=triangnew(primo,2)
                 ENDIF
                 IF (sqrt((xRnew(nodo4)-xcirccentrnew(primo))**2     &
                      &  +(yRnew(nodo4)-ycirccentrnew(primo))**2).LT. &
                      &  raggiocirccentrnew(primo)*(1.d0-eps1)) THEN
                    triangnew4(1)=nodo2
                    triangnew4(2)=seg(primo,j,1)
                    triangnew4(3)=nodo4
                    CALL nuovotriangolo2(triangnew4,segnew4, &
                         &  xRnew(nodo2),yRnew(nodo2),          &
                         &  xRnew(seg(primo,j,1)),yRnew(seg(primo,j,1)),&
                         &  xRnew(nodo4),yRnew(nodo4),xcirco,ycirco, &
                         &  raggiocirco)
                    triangnew(primo,1)=triangnew4(1)
                    triangnew(primo,2)=triangnew4(2)
                    triangnew(primo,3)=triangnew4(3)
                    triangnew4(1)=nodo4
                    triangnew4(2)=nodo2
                    triangnew4(3)=seg(primo,j,2)
                    xcirccentrnew(primo)=xcirco
                    ycirccentrnew(primo)=ycirco
                    raggiocirccentrnew(primo)=raggiocirco
                    CALL nuovotriangolo2(triangnew4,segnew5, &
                         &  xRnew(nodo4),yRnew(nodo4),          &
                         &  xRnew(nodo2),yRnew(nodo2),          &
                         &  xRnew(seg(primo,j,2)),yRnew(seg(primo,j,2)),&
                         &  xcirco,ycirco,raggiocirco)
                    triangnew(cambia(i),1)=triangnew4(1)
                    triangnew(cambia(i),2)=triangnew4(2)
                    triangnew(cambia(i),3)=triangnew4(3)
                    xcirccentrnew(cambia(i))=xcirco
                    ycirccentrnew(cambia(i))=ycirco
                    raggiocirccentrnew(cambia(i))=raggiocirco
                    DO l=1,3
                       seg(primo,l,1)=segnew4(l,1)
                       seg(primo,l,2)=segnew4(l,2)
                    ENDDO
                    DO l=1,3
                       seg(cambia(i),l,1)=segnew5(l,1)
                       seg(cambia(i),l,2)=segnew5(l,2)
                    ENDDO
                    ntg=0
                    cambia(quanti+1)=cambia(1)
                    cambia(quanti+2)=cambia(i)
                    DO l=1,i-2
                       cambia(l)=cambia(l+1)
                    ENDDO
                    DO l=i-1,quanti
                       cambia(l)=cambia(l+2)
                    ENDDO
                    GOTO 800
                 ENDIF
              ENDIF
           ENDDO
        ENDDO
     ENDDO
     cambia(quanti+1)=cambia(1)
     DO l=1,quanti
        cambia(l)=cambia(l+1)
     ENDDO
     ntg=ntg+1
     GOTO 800
  ENDIF
  
313 CONTINUE
  
  DO i=1,quanti
     triangnew4(1)=triangnew(cambia(i),1)
     triangnew4(2)=triangnew(cambia(i),2)
     triangnew4(3)=triangnew(cambia(i),3)
     CALL nuovotriangolo(triangnew4,segnew5,xRnew(triangnew4(1)), &
          &  yRnew(triangnew4(1)),xRnew(triangnew4(2)),yRnew(triangnew4(2)) &
          &  ,xRnew(triangnew4(3)),yRnew(triangnew4(3)),   &
          &  xcirco,ycirco,raggiocirco,                    &
          &  xbaric(cambia(i)),ybaric(cambia(i)))
     
     IF ((triangnew4(1).EQ.2).AND.(triangnew4(2).EQ.3).AND. &
          & (triangnew4(3).EQ.20)) THEN
        
     ENDIF
     
  ENDDO
  
  RETURN
END SUBROUTINE swappalista

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!                       !!!!!!!!!!!
!!!!!!!!!   CERCATRIANGOLIVICINI     !!!!!!
!!!!!!!!!                       !!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE cercatriangolivicini(npox,ntrix,triang,nt,quale,quantivicini,&
     &  vicini,seg,nodoopposto)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: npox,ntrix ! max number of points and triangles
  INTEGER, INTENT(IN) :: triang(ntrix,4)
  INTEGER, INTENT(IN) :: nt,quale,seg(ntrix,3,2)
  INTEGER, INTENT(OUT) :: vicini(3),nodoopposto(3),quantivicini
! END INTERFACE
! =========================================================================
  INTEGER h,i,j
  
  quantivicini=0
  DO 1000 i=1,nt
     DO j=1,3
        DO h=1,3
           IF ((seg(quale,j,1).EQ.seg(i,h,1)).AND.               &
                &(seg(quale,j,2).EQ.seg(i,h,2)).AND.(i.ne.quale))THEN
              quantivicini=quantivicini+1
              IF (h.EQ.1) THEN
                 nodoopposto(quantivicini)=triang(i,3)
              ELSEIF (h.EQ.2) THEN
                 nodoopposto(quantivicini)=triang(i,1)
              ELSE
                 nodoopposto(quantivicini)=triang(i,2)
              ENDIF
              vicini(quantivicini)=i
              GOTO 1000
           ENDIF
        ENDDO
     ENDDO
1000 ENDDO
  RETURN
END SUBROUTINE cercatriangolivicini

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!                       !!!!!!!!!!!
!!!!!!!!!       VERIFICAINT         !!!!!!!
!!!!!!!!!                          !!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE verificaint(npox,ntrix,triang,nt,xR5,yR5,intersez2,&
     & xcentro,ycentro,raggio,nodibad,noditot,triangnew,segmento)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: npox,ntrix ! max number of points and triangles
  INTEGER, INTENT(IN) :: nt,noditot
  INTEGER, INTENT(OUT) :: intersez2,nodibad
  INTEGER, INTENT(IN) :: triangnew(4),triang(ntrix,4)
  DOUBLE PRECISION, INTENT(IN) :: xR5(0:npox),yR5(0:npox)
  INTEGER, INTENT(IN) :: segmento(npox,npox)
  ! END INTERFACE
  DOUBLE PRECISION :: x1,x2,x3,x4,ba11,ba12,ba21,ba22,bb1,bb2
  DOUBLE PRECISION :: y1,y2,y3,y4,lambda1,gamma1
  INTEGER :: i,j,h ! loop indexes
  DOUBLE PRECISION :: xcentro,ycentro,raggio
! =========================================================================
  
  intersez2=0
  DO i=1,noditot-1
     DO j=i,noditot
        IF ((segmento(i,j).ne.0).AND.&
             &(((i.ne.triangnew(1)).OR.(i.ne.triangnew(3))).OR.&
             &((j.ne.triangnew(3)).OR.(j.ne.triangnew(1))))) THEN
           x1=xR5(triangnew(1))
           x2=xR5(triangnew(3))
           x3=xR5(i)
           x4=xR5(j)
           y1=yR5(triangnew(1))
           y2=yR5(triangnew(3))
           y3=yR5(i)
           y4=yR5(j)
           ba11=x1-x2
           ba12=x4-x3
           ba21=y1-y2
           ba22=y4-y3
           IF (abs(ba11*ba22-ba12*ba21)/(sqrt(ba11**2+ba21**2)&
                & *sqrt(ba12**2+ba22**2)).ge.1.d-4) THEN
              bb1=x4-x2
              bb2=y4-y2
              lambda1=(bb1*ba22-bb2*ba12)/(ba11*ba22-ba12*ba21)
              gamma1=(ba11*bb2-ba21*bb1)/(ba11*ba22-ba12*ba21)
           ELSE
              lambda1=2
           ENDIF
           IF (((lambda1.GT.1.d-4).AND.(lambda1.LT.1-1.d-4)).AND.&
                & ((gamma1.GT.1.d-4).AND.(gamma1.LT.1-1.d-4))) THEN
              
              intersez2=1
              GOTO 2000
           ENDIF
        ENDIF
        IF ((segmento(i,j).ne.0).AND.&
             &(((i.ne.triangnew(2)).OR.(i.ne.triangnew(3))).OR.&
             &((j.ne.triangnew(3)).OR.(j.ne.triangnew(2))))) THEN
           x1=xR5(triangnew(2))
           x2=xR5(triangnew(3))
           x3=xR5(i)
           x4=xR5(j)
           y1=yR5(triangnew(2))
           y2=yR5(triangnew(3))
           y3=yR5(i)
           y4=yR5(j)
           ba11=x1-x2
           ba12=x4-x3
           ba21=y1-y2
           ba22=y4-y3
           IF (abs(ba11*ba22-ba12*ba21)/(sqrt(ba11**2+ba21**2) &
                & *sqrt(ba12**2+ba22**2)).ge.1.d-4) THEN
              bb1=x4-x2
              bb2=y4-y2
              lambda1=(bb1*ba22-bb2*ba12)/(ba11*ba22-ba12*ba21)
              gamma1=(ba11*bb2-ba21*bb1)/(ba11*ba22-ba12*ba21)
           ELSE
              lambda1=2
           ENDIF
           IF (((lambda1.GT.1.d-4).AND.(lambda1.LT.1-1.d-4)).AND.&
                &((gamma1.GT.1.d-4).AND.(gamma1.LT.1-1.d-4))) THEN
              intersez2=1
              
              GOTO 2000
           ENDIF
        ENDIF
     ENDDO
  ENDDO
  IF (nt.ge.2) THEN
     DO i=1,nt-1
        IF (((triang(i,1).EQ.triangnew(1)).AND.&
             &(triang(i,2).EQ.triangnew(2))).AND.&
             &(triang(i,3).EQ.triangnew(3))) THEN
           intersez2=1
           GOTO 2000
        ENDIF
     ENDDO
     
  ENDIF
  
  DO i =1,noditot
     IF((i.ne.triangnew(1)).AND.(i.ne.triangnew(2)).AND.&
          &(i.ne.triangnew(3))) THEN
        x1 = xR5(triangnew(1))
        x2 = xR5(triangnew(2))
        x3 = xR5(triangnew(3))
        x4 = xR5(i)
        y1 = yR5(triangnew(1))
        y2 = yR5(triangnew(2))
        y3 = yR5(triangnew(3))
        y4 = yR5(i)
        
        ba11=x2-x1
        ba12=x3-x1
        ba21=y2-y1
        ba22=y3-y1
        IF (abs(ba11*ba22-ba12*ba21)/(sqrt(ba11**2+ba21**2)&
             & *sqrt(ba12**2+ba22**2)).ge.1.d-4) THEN
           bb1=x4-x1
           bb2=y4-y1
           lambda1=(bb1*ba22-bb2*ba12)/(ba11*ba22-ba12*ba21)
           gamma1=(ba11*bb2-ba21*bb1)/(ba11*ba22-ba12*ba21)
        ELSE
           intersez2=1
           GOTO 2000
        ENDIF
        
        IF (((lambda1.GT.-1.d-4).AND.(lambda1.LT.1+1.d-4)).AND.&
             & ((gamma1.GT.-1.d-4).AND.(gamma1.LT.1+1.d-4)).AND.&
             &(lambda1+gamma1.le.1.E0+1.d-4)) THEN
           intersez2=1
           
           GOTO 2000
        ENDIF
     ENDIF
  ENDDO
  
  nodibad=0
  DO i=1,noditot
     IF (((i.EQ.triangnew(1)).OR.(i.EQ.triangnew(2)))&
          &.OR.(i.EQ.triangnew(3))) THEN
        nodibad=nodibad+1
     ELSE
        IF ((xcentro-xR5(i))**2+(ycentro-yR5(i))**2.le.raggio**2) &
             &THEN
           nodibad=nodibad+1
        ENDIF
     ENDIF
  ENDDO
2000 CONTINUE
  
  RETURN     
END SUBROUTINE verificaint

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!			     !!!!!!!!!!!
!!!!!!!!!   CERCANODIVICINI          !!!!!!!!!!!
!!!!!!!!!			     !!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE cercanodivicini(npox,ntrix,quale,seg,xR,yR,density,nt)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: npox,ntrix ! max number of points and triangles
  INTEGER,INTENT(IN) :: nt
  INTEGER,INTENT(IN) :: quale,seg(2*npox,3,2)
  DOUBLE PRECISION,INTENT(INOUT) :: xR(0:npox),yR(0:npox)
  DOUBLE PRECISION,INTENT(INOUT) :: density(npox)
  ! end interface
  DOUBLE PRECISION :: xnew,ynew,sum
  DOUBLE PRECISION :: densityNew
  INTEGER :: i,j,h
  INTEGER :: vertice(npox),newlist(npox),vertmax
  INTEGER :: indmax
! ==========================================================
  
  h=0
  DO i=1,nt
     DO j=1,3
        IF (seg(i,j,1).EQ.quale) THEN
           h=h+1
           vertice(h)=seg(i,j,2)
        ELSEIF (seg(i,j,2).EQ.quale) THEN
           h=h+1
           vertice(h)=seg(i,j,1)
        ENDIF
     ENDDO
  ENDDO
  
  DO i=1,h
     vertmax=0
     DO j=i,h
        IF (vertice(j).GT.vertmax) THEN
           vertmax=vertice(j)
           indmax=j
        ENDIF
     ENDDO
     vertice(indmax)=vertice(i)
     vertice(i)=vertmax
  ENDDO
  newlist(1)=vertice(1)
  j=1
  DO i=2,h
     IF (vertice(i).ne.newlist(j)) THEN
        j=j+1
        newlist(j)=vertice(i)
     ENDIF
  ENDDO
  xnew=0.
  ynew=0.
  sum=0.
  densityNew=0.
  DO i=1,j
     xnew=xnew+1/density(newlist(i))*xR(newlist(i))
     ynew=ynew+1/density(newlist(i))*yR(newlist(i))
     sum=sum+1/density(newlist(i))
     
  ENDDO
  xnew=xnew/sum
  ynew=ynew/sum
  densityNew=j/sum
  xR(quale)=xnew
  yR(quale)=ynew
  density(quale)=densityNew
END SUBROUTINE cercanodivicini

! *************************************************
! SMOOTHING and OPTIMIZATION of the TRIANGULATION
! *************************************************
SUBROUTINE smoothing(npox,ntrix,eps,nodiest,noditot,nt,seg, &
     & xR,yR,density,triang)
  IMPLICIT none
  INTEGER, INTENT(IN) :: npox,ntrix ! max number of points and triangles
  DOUBLE PRECISION,INTENT(IN) :: eps
  INTEGER,INTENT(IN) :: nodiest,noditot,nt
  INTEGER,INTENT(INOUT) :: seg(ntrix,3,2)
  DOUBLE PRECISION,INTENT(INOUT) :: xR(0:npox),yR(0:npox)
  DOUBLE PRECISION,INTENT(INOUT) :: density(npox)
  INTEGER,INTENT(INOUT) :: triang(ntrix,4)
! END INTERFACE
  INTEGER :: i,j,l,h ! loop indexes     
  INTEGER :: count ! counter
  INTEGER :: triangnew(4)
  INTEGER :: segnew(3,2) !only for CALL od nuovotriangolo, not needed here
  DOUBLE PRECISION :: xcirccentr(2*npox),ycirccentr(2*npox), &
       & raggiocirccentr(2*npox)
  DOUBLE PRECISION :: xcentro,ycentro,raggiont
  INTEGER :: nodoopposto(3),quantivicini,vicini(3),quanti2,cambia(2*npox)
  DOUBLE PRECISION :: xbaric(2*npox),ybaric(2*npox)
! ================================================================
  
! move the nodes in the barycenter of the neighbor polygon
  DO h=1,2
     DO 350 l=nodiest+1,noditot
        CALL cercanodivicini(npox,ntrix,l,seg,xR,yR,density,nt)
350  ENDDO
     DO 360 l=noditot,nodiest+1,-1
        CALL cercanodivicini(npox,ntrix,l,seg,xR,yR,density,nt)
360  ENDDO
  ENDDO
  
  DO i=1,nt
     
     triangnew(1)=triang(i,1)
     triangnew(2)=triang(i,2)
     triangnew(3)=triang(i,3)
     triangnew(4)=triang(i,1)
     CALL nuovotriangolo(triangnew,segnew, &
          &  xR(triangnew(1)),yR(triangnew(1)), &
          &  xR(triangnew(2)),yR(triangnew(2)), &
          &  xR(triangnew(3)),yR(triangnew(3)), &
          &  xcentro,ycentro,raggiont, &
          &  xbaric(i),ybaric(i))
     xcirccentr(i)=xcentro
     ycirccentr(i)=ycentro
     raggiocirccentr(i)=raggiont
     
  ENDDO
  
  count = 0
  
400 count = count+1
  !     if (count.GT.2**nt) THEN
  IF(count.gt.2*nt)THEN
     GOTO 133
  ENDIF
  
  DO i=1,nt
     CALL cercatriangolivicini(npox,ntrix,triang,nt,i,quantivicini, &
          &  vicini,seg,nodoopposto)
     
     DO j=1,quantivicini
        IF (sqrt((xR(nodoopposto(j))-xcirccentr(nt))**2  &
             &  +(yR(nodoopposto(j))-ycirccentr(nt))**2)   &
             &  .LT.raggiocirccentr(nt)*(1.d0-eps)) THEN
           CALL listaswap(npox,ntrix,triang,xR,yR,xcirccentr, &
                &   seg,ycirccentr,raggiocirccentr, &
                &   nt,quanti2,cambia,nt,eps)
           CALL swappalista(npox,ntrix,triang,xR,yR,xcirccentr,seg, &
                &   ycirccentr,raggiocirccentr,quanti2,cambia, &
                &   xbaric,ybaric,eps)
           
           GOTO 400
           
        ENDIF
     ENDDO
  ENDDO
  
133 CONTINUE 
  
END SUBROUTINE smoothing

! =============================
! READING TRIANGULATION FILES
! =========================================================================
SUBROUTINE outtriang(trifil,noditot,nt,xR,yR,triang)
  CHARACTER*(*),INTENT(IN) :: trifil
  INTEGER, INTENT(IN) :: noditot,nt,triang(nt,4)
  DOUBLE PRECISION, INTENT(IN) :: xR(noditot),yR(noditot)
  ! END INTERFACE
! =========================================================================
  DOUBLE PRECISION :: xn1,xn2,xn3,yn1,yn2,yn3
  INTEGER :: i
  INTEGER :: iun ! unit number      
  CALL filopn(iun,trifil,'unknown')  
  !       write total number of nodes and of triangles
  write(iun,100) noditot,nt,nfunc
100 FORMAT(i5,8x,i5,8x,i5)
  DO i=1,noditot
     write(iun,101) xR(i),yR(i), 0.d0
  ENDDO
101 FORMAT(f17.12,3x,f17.12,3x,f3.1)
  DO i=1,nt
     xn1=xR(triang(i,1))
     yn1=yR(triang(i,1))
     xn2=xR(triang(i,2))
     yn2=yR(triang(i,2))
     xn3=xR(triang(i,3))
     yn3=yR(triang(i,3))    
     IF ((xn2-xn1)*(yn3-yn1)-(xn3-xn1)*(yn2-yn1).LT.0) THEN
        write(iun,100) triang(i,1),triang(i,3),triang(i,2)
     ELSE
        write(iun,100) triang(i,1),triang(i,2),triang(i,3)
     ENDIF
  ENDDO
  CALL filclo(iun,' ')
END SUBROUTINE outtriang

SUBROUTINE intriang(npox,ntrix,trifil,noditot,nt,xR,yR,triang)
!  USE triangles
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: npox,ntrix ! max number of points and triangles
  CHARACTER*(*),INTENT(IN) :: trifil ! file name
  INTEGER, INTENT(OUT) :: noditot,nt ! total number of nodes and triangles
  INTEGER, INTENT(OUT) :: triang(3,ntrix) ! incidence relations between nodes
  DOUBLE PRECISION, INTENT(OUT) :: xR(npox),yR(npox) ! node coords
  ! END INTERFACE
! =========================================================================
  INTEGER :: i ! loop index
  INTEGER :: iun ! unit number      
  CALL filopn(iun,trifil,'old')  
  ! READ total number of nodes and of triangles
  READ(iun,100) noditot,nt
100 FORMAT(i5,8x,i5,8x,i5)
  ! READ  node coordinates
  DO i=1,noditot
     READ(iun,101) xR(i),yR(i)
  enddo
101 FORMAT(f17.12,3x,f17.12,3x,f3.1)
  ! READ incidence relations between nodes
  DO  i=1,nt
     READ(iun,100) triang(1,i),triang(2,i),triang(3,i)
  ENDDO
  CALL filclo(iun,' ')
END SUBROUTINE intriang

! =================================================================
! MERGE_TRIANG
! merges two triangulations into one, with list of nodes and list
! of triangles
! =================================================================
SUBROUTINE merge_triang(npox,ntrix,frv1,rdot1,npotot1,triang1,ntri1, &
     & frv2,rdot2,npotot2,triang2,ntri2, &
     & logrv,rdot,npotot,triang,ntri)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: npox,ntrix ! max number of points and triangles
! triangulations: node coordinates (log10(r), rdot)
  DOUBLE PRECISION,INTENT(IN),DIMENSION(npox) :: frv1,frv2,rdot1,rdot2
  DOUBLE PRECISION,INTENT(OUT),DIMENSION(npox) :: logrv,rdot
! triangle indexes
  INTEGER,INTENT(IN) :: triang1(ntrix,4),triang2(ntrix,4)
  INTEGER,INTENT(OUT) :: triang(ntrix,4)
! total number of points
  INTEGER,INTENT(IN) :: npotot1,npotot2 
  INTEGER,INTENT(OUT) :: npotot
! total number of triangles
  INTEGER,INTENT(IN) :: ntri1,ntri2
  INTEGER,INTENT(OUT) :: ntri
! =========== end interface =======================================
! to count the node repetitions
  INTEGER nrep,repnod2list(npotot2),repnod1list(npotot2)
! loop indexes
  INTEGER i,j,h,m,n
! =================================================================

! initialization
  nrep=0
  
! nodes and triangles of the first list remain unchanged
  DO i = 1,npotot1
     logrv(i) = frv1(i)
     rdot(i) = rdot1(i)
     
  ENDDO
  DO i=1,ntri1
     DO n = 1,3
        triang(i,n) = triang1(i,n)
     ENDDO
  ENDDO
  
  DO h = 1,npotot2
     DO i = 1,npotot1
        
        IF (abs((frv2(h) - frv1(i))/frv2(h)).lt.0.d0) THEN
           IF (abs((rdot2(h) - rdot1(i))).lt.0.d0) THEN
                  ! updating counters
              nrep = nrep + 1       ! count node repetitions
! number of the repeated node in 2nd list 
              repnod2list(nrep) = h
! number of the repeated node in 1st list
              repnod1list(nrep) = i
!                  WRITE(*,*)' doppio ', nrep,h,i                  
              GOTO 10
           ELSE
! add to the list as a new node
              logrv(h+npotot1-nrep) = frv2(h)
              rdot(h+npotot1-nrep) = rdot2(h)
           ENDIF
        ELSE
! add to the list as a new node
           logrv(h+npotot1-nrep) = frv2(h)
           rdot(h+npotot1-nrep) = rdot2(h)
        ENDIF
            
     ENDDO
         
10   CONTINUE
  ENDDO

! LOOP ON TRIANG      
  DO i = 1,ntri2
     DO n = 1,3
! no duplicates
        IF(nrep.eq.0)THEN
           triang(i+ntri1,n) = triang2(i,n)+npotot1
        ELSE
! there are duplicates            
           IF (triang2(i,n).lt.repnod2list(1)) THEN
              triang(i+ntri1,n) = triang2(i,n)+npotot1
           ENDIF
           
           DO m =1,nrep-1
              IF (triang2(i,n).eq.repnod2list(m)) THEN
                 triang(i+ntri1,n) = repnod1list(m)
              ELSEIF ((triang2(i,n).gt.repnod2list(m)).and. &
                   &             (triang2(i,n).lt.repnod2list(m+1))) THEN
                 triang(i+ntri1,n) = triang2(i,n)+npotot1-m
              ENDIF
           ENDDO
            
           IF (triang2(i,n).eq.repnod2list(nrep)) THEN
              triang(i+ntri1,n) = repnod1list(nrep)
           ELSEIF (triang2(i,n).gt.repnod2list(nrep)) THEN
              triang(i+ntri1,n) = triang2(i,n)+npotot1-nrep
           ENDIF
        ENDIF
     ENDDO
  ENDDO

! total number of triangles
  ntri = ntri1+ntri2
  npotot = npotot1 + npotot2 - nrep      

! writing on output file
  OPEN(unit=2,file='merge_tri.dat',status='unknown') 
  WRITE(2,100) npotot,ntri,0
  DO i=1,npotot 
     WRITE(2,101) logrv(i),rdot(i),0.d0
  ENDDO
  DO j=1,ntri
     WRITE(2,100) triang(j,1),triang(j,2),triang(j,3)
  ENDDO
  
100 FORMAT(i5,8x,i5,8x,i5)
101 FORMAT(f10.5,3x,f10.5,3x,f10.5)
  
  CLOSE(2)
  
  RETURN
END SUBROUTINE merge_triang

! ================================================
!   written by MATTIA DE MICHIELI VITTURI 2002
!       =======================================
!       =======================================
! last modified GFG 7/7/2003
SUBROUTINE uniform(xRR,yRR,nodiinit,nodifin,quali,npoxmul)
  implicit none
  INTEGER, INTENT(IN) :: nodiinit,nodifin,npoxmul
  INTEGER, INTENT(OUT) :: quali(npoxmul) 
  DOUBLE PRECISION ,INTENT(IN) :: xRR(nodiinit),yRR(nodiinit)      
!     END INTERFACE
  DOUBLE PRECISION :: lung(0:nodiinit),densita(nodiinit),lung2
  DOUBLE PRECISION :: lungtot,distmin,distanza(nodiinit)
  INTEGER :: nodoelim
  INTEGER :: quanti,i,j,h,dist
  INTEGER :: ninit,nfin
  DOUBLE PRECISION :: intervallo
  DOUBLE PRECISION, PARAMETER ::opt=1.d0
! ================================================      
  ninit = nodiinit
  nfin = nodifin
  quanti=ninit
  
  DO i=1,quanti
     quali(i)=i
  ENDDO
      
!      IF (ninit.eq.2) THEN
!         write(*,*)'xRR,yRR',xRR,yRR
!         write(*,*)'ninit,nfin',ninit,nfin
!      ENDIF
  IF(ninit.lt.nfin) THEN
     write(*,*)'ninit,nfin',ninit,nfin
     write(*,*)'ninit < nfin: forcing ninit = nfin'
     nfin = ninit
  ENDIF
  IF(nfin.lt.2) THEN
     IF(ninit.ge.2) THEN
        write(*,*)'nfin=',nfin,'  forcing nfin = 2'
        nfin = 2
     ELSE
        WRITE(*,*)'ERROR: ninit,nfin',ninit,nfin
     ENDIF
  ENDIF
! ======================================
! =======  CALCOLA ASC.CURV.  ==========
! ======================================
  lung(0)=0.d0
  lungtot=0.d0
  DO i = 1,ninit-1
     lung(i)=sqrt((xRR(i)-xRR(i+1))**2+(yRR(i)-yRR(i+1))**2)
     lungtot=lungtot+lung(i)
  ENDDO
  lung(ninit)=0.
  intervallo=(lungtot/(nfin-1))  
! ======================================
! ======  CALCOLA DIST E DENS  =========
! ======================================
  densita(1)=0.
  distanza(1)=0.
  DO i = 2,ninit-1
     densita(i)=opt*min(lung(i),lung(i-1))+(1-opt)*(lung(i)+lung(i-1))
     distanza(i)=lungtot
     lung2=0.
     DO j= 1,i
        lung2=lung2+lung(j-1)
     ENDDO
     DO j= 1,nfin-2
        distanza(i)=min(distanza(i),abs(lung2-lungtot*j/(nfin-1)))
     ENDDO
  ENDDO
  densita(ninit)=0.
  distanza(ninit)=0.  
! ======================================
! =========  ELIMINA I PUNTI  ==========
! ======================================
  DO h=1,ninit-nfin
     distmin=2*lungtot
     DO i=2,quanti-1
        IF (densita(i)/(1+sqrt(distanza(i))).lt.distmin) THEN
           distmin=densita(i)/(1+sqrt(distanza(i)))
           nodoelim=i   !indice punto da eliminare
        ENDIF
     ENDDO
     IF (nodoelim.eq.2) THEN
        lung(1)=lung(1)+lung(2)
        densita(3)=opt*min(lung(1),lung(3))+&
             &  (1-opt)*(lung(1)+lung(3))
        DO j=nodoelim,quanti-2
           quali(j)=quali(j+1)
           densita(j)=densita(j+1)
           lung(j)=lung(j+1)
           distanza(j)=distanza(j+1) 
        ENDDO
     ELSEIF (nodoelim.eq.quanti-1) THEN
        lung(quanti-2)=lung(quanti-2)+lung(quanti-1)
        densita(quanti-2)=opt*(min(lung(quanti-3),lung(quanti-2)))+&
             &  (1-opt)*(lung(quanti-3)+lung(quanti-2))
        quali(quanti-1)=quali(quanti)
     ELSE
        lung(nodoelim-1)=lung(nodoelim-1)+lung(nodoelim)
        densita(nodoelim-1)=opt*min(lung(nodoelim-2),&
             &   lung(nodoelim-1))+&
             &   (1-opt)*(lung(nodoelim-2)+lung(nodoelim-1))
        DO j=nodoelim,quanti-2
           quali(j)=quali(j+1)
           lung(j)=lung(j+1)
           densita(j)=opt*min(lung(j),lung(j-1))+&
                &   (1-opt)*(lung(j)+lung(j-1))
           distanza(j)=distanza(j+1) !!
        ENDDO
     ENDIF
     quali(quanti-1)=quali(quanti)
     quanti=quanti-1
  ENDDO
  
END SUBROUTINE uniform
