i7trivox1.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    i7trivox1                  ../starter/source/interfaces/inter3d1/i7trivox1.F
!||--- called by ------------------------------------------------------
!||    i7buc_vox1                 ../starter/source/interfaces/inter3d1/i7buc_vox1.F
!||--- calls      -----------------------------------------------------
!||    i7cor3                     ../starter/source/interfaces/inter3d1/i7cor3.F
!||    i7dst3                     ../starter/source/interfaces/inter3d1/i7dst3.F
!||    i7pen3                     ../starter/source/interfaces/inter3d1/i7pen3.F
!||    inter_save_candidate       ../starter/source/interfaces/inter3d1/inter_save_candidate.F90
!||--- uses       -----------------------------------------------------
!||    inter_save_candidate_mod   ../starter/source/interfaces/inter3d1/inter_save_candidate.F90
!||    tri7box                    ../starter/share/modules1/tri7box.F
!||====================================================================
      SUBROUTINE i7trivox1(
     1      NSN    ,I_MEM   ,IRECT    ,X       ,STF     ,
     2      STFN   ,XYZM    ,NSV      ,
     3      MULNSN ,NOINT   ,TZINF    ,GAP_S_L ,GAP_M_L ,
     4      VOXEL  ,NBX     ,NBY      ,NBZ     ,NRTM_L    ,
     5      IGAP   ,GAP     ,GAP_S    ,GAP_M   ,GAPMIN  ,
     6      GAPMAX ,MARGE   ,CURV_MAX ,BGAPSMX ,ISTF    ,
     7      I_STOK ,nin,
     8      ID     ,TITR    ,DRAD     ,INDEX   ,
     9      IREMNODE,FLAGREMNODE,KREMNODE,REMNODE,
     1      DGAPLOAD,ipari,intbuf_tab,
     2      iix,iiy,iiz,local_next_nod,nrtm,IS_USED_WITH_LAW151 )
C============================================================================
C   M o d u l e s
C-----------------------------------------------
      USE my_alloc_mod
      USE tri7box
      USE names_and_titles_mod , ONLY : nchartitle
      use inter_save_candidate_mod , only : inter_save_candidate
      use array_mod 
      use intbufdef_mod 
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
c     parameter setting the size for the vector (orig version is 128)
      INTEGER NVECSZ
      parameter(nvecsz = mvsiz)
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com04_c.inc"
#include      "param_c.inc"
C-----------------------------------------------
C   ROLE OF THE ROUTINE:
C   ===================
C   CLASSIFIES NODES INTO BOXES
C   SEARCHES FOR EACH FACETTE OF CONCERNED BOXES
C   SEARCH FOR CANDIDATES
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C
C     NOM          DESCRIPTION                       E/S
C
C     IRECT(4,*)   ARRAY OF CONEC FACETTES        E
C     X(3,*)       COORDONNEES NODALES               E
C     NSV          NOS SYSTEMES DES NODES E
C     Xmax larger abcisse existing e
C     XMAX largest order.existing E
C     Xmax larger existing side E
C     I_STOK storage level of couples
C                                CANDIDATES impact    E/S
C     CAND_N       boites resultats nodes C     CAND_E       adresses des boites resultat elements
C                  MULNSN = MULTIMP*NSN MAX SIZE NOW ALLOWED FOR
C                  COUPLES NODES,ELT CANDIDATES
C     NOINT        INTERFACE USER NUMBER
C     TZINF        TAILLE ZONE INFLUENCE
C
C     Prov_n Provisional Cand_n (static variable in i7tri)
C     PROV_E       CAND_E provisoire (variable static in i7tri)
 
C     VOXEL(ix,iy,iz) contains the local number of the first node of
C                  the box
C     Local_next_nod (I) Next node in the same box (if /= 0)
C     Last_nod (i) Last node in the same box (if /= 0)
C                  used only to go directly from the first
C                       node at the last
C     INDEX        index of IRECTM active for this domain
C     NRTM_L         number of active segments for this domain
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      integer, intent(in) :: nrtm !< number of segment
      integer, intent(in) :: nrtm_l !< number of segment on the current proc
      INTEGER I_MEM,NSN,
     .        MULNSN,NOINT,IGAP,NBX,NBY,NBZ,IREMNODE,FLAGREMNODE,
     .        NSV(NSN),
     .        IRECT(4,NRTM), VOXEL(NBX+2,NBY+2,NBZ+2),ISTF,
     .        I_STOK ,J_STOK,
     .        INDEX(*),KREMNODE(*),REMNODE(*)
      my_real
     .   X(3,*),XYZM(6,2),STF(*),STFN(*),GAP_S(*),GAP_M(*),
     .   TZINF,MARGE,GAP,GAPMIN,GAPMAX,BGAPSMX,DRAD,
     .   curv_max(*),gap_s_l(*),gap_m_l(*)
      INTEGER ID
      CHARACTER(LEN=NCHARTITLE)::TITR
      integer, intent(in) :: nin !< interface index
      integer, dimension(npari), intent(inout) :: ipari !< interface data
      type(intbuf_struct_), intent(inout) :: intbuf_tab !< interface data
      INTEGER,  dimension(nsn), intent(inout) :: iix,iiy,iiz,local_next_nod
      LOGICAL,INTENT(IN) :: IS_USED_WITH_LAW151
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER NB_NCN,NB_NCN1,NB_ECN,I,J,DIR,NB_NC,NB_EC,
     .        n1,n2,n3,n4,nn,ne,k,l,ncand_prov,ii,jj,kk,
     .        nsnf, nsnl, i_bid,delnod
      my_real
     .   dx,dy,dz,xs,ys,zs,xx,sx,sy,sz,s2,
     .   xmin, xmax,ymin, ymax,zmin, zmax, tz, gapsmx, gapl,
     .   xx1,xx2,xx3,xx4,yy1,yy2,yy3,yy4,zz1,zz2,zz3,zz4,
     .   d1x,d1y,d1z,d2x,d2y,d2z,dd1,dd2,d2,a2,gs,gapv(mvsiz)
      INTEGER  LAST_NOD(NSN)
      INTEGER*8 meanjj8
      INTEGER  IX,IY,IZ,NEXT,M1,M2,M3,M4,
     .         ix1,iy1,iz1,ix2,iy2,iz2
      integer, dimension(mvsiz) :: prov_n,prov_e
      my_real
     .   xminb,yminb,zminb,xmaxb,ymaxb,zmaxb,
     .   xmine,ymine,zmine,xmaxe,ymaxe,zmaxe,aaa,tstart,tstop
      my_real , INTENT(IN) :: dgapload
      INTEGER FIRST,NEW,LAST,M
      INTEGER, DIMENSION(MVSIZ) :: IX11,IX12,IX13,IX14,NSVG
      my_real, DIMENSION(MVSIZ) :: x1,x2,x3,x4
      my_real, DIMENSION(MVSIZ) :: y1,y2,y3,y4
      my_real, DIMENSION(MVSIZ) :: z1,z2,z3,z4
      my_real, DIMENSION(MVSIZ) :: xi,yi,zi
      my_real, DIMENSION(MVSIZ) :: x0,y0,z0
      my_real, DIMENSION(MVSIZ) :: nx1,ny1,nz1
      my_real, DIMENSION(MVSIZ) :: nx2,ny2,nz2
      my_real, DIMENSION(MVSIZ) :: nx3,ny3,nz3
      my_real, DIMENSION(MVSIZ) :: nx4,ny4,nz4
      my_real, DIMENSION(MVSIZ) :: p1,p2,p3,p4
      my_real, DIMENSION(MVSIZ) :: lb1,lb2,lb3,lb4
      my_real, DIMENSION(MVSIZ) :: lc1,lc2,lc3,lc4
      my_real, DIMENSION(MVSIZ) :: n11,n21,n31
      my_real, DIMENSION(MVSIZ) :: stif,pene
 
      INTEGER, DIMENSION(:), ALLOCATABLE ::  TAGNOD
 
      integer , external :: omp_get_thread_num,omp_get_num_threads
      integer :: itask,nthreads
      integer :: my_old_size,my_address
      integer :: local_i_stok,multimp
      integer :: local_cand_n_size,local_cand_e_size
 
      integer, save :: i_stok_old
      integer, dimension(:), allocatable, save :: cand_n_size,cand_e_size
      integer, dimension(:), allocatable, save :: address_cand_n,address_cand_e
      type(array_type_int_1d) :: local_cand_n
      type(array_type_int_1d) :: local_cand_e
 
      integer :: my_size,mode
      integer, dimension(:), allocatable :: my_index
      integer, dimension(:,:), allocatable :: sort_array,save_array
      integer, dimension(70000) :: work
      integer stagnod
C-----------------------------------------------
      ! allocation of local arrays 
      itask = omp_get_thread_num()
      nthreads = omp_get_num_threads()
      local_cand_n_size = size(intbuf_tab%cand_n) / nthreads
      local_cand_e_size = size(intbuf_tab%cand_e) / nthreads
      local_i_stok = 0
      local_cand_n%size_int_array_1d = local_cand_n_size
      local_cand_e%size_int_array_1d = local_cand_e_size
      call alloc_1d_array(local_cand_n)
      call alloc_1d_array(local_cand_e)
 
      xmin = xyzm(1,1)
      ymin = xyzm(2,1)
      zmin = xyzm(3,1)
      xmax = xyzm(4,1)
      ymax = xyzm(5,1)
      zmax = xyzm(6,1)
 
      xminb = xyzm(1,2)
      yminb = xyzm(2,2)
      zminb = xyzm(3,2)
      xmaxb = xyzm(4,2)
      ymaxb = xyzm(5,2)
      zmaxb = xyzm(6,2)
 
C=======================================================================
C 1 placing nodes in boxes
C=======================================================================
 
!$OMP MASTER
 
      allocate( cand_n_size(nthreads+1),cand_e_size(nthreads+1) )
      allocate( address_cand_n(nthreads+1),address_cand_e(nthreads+1) )
      cand_n_size(1:nthreads+1) = 0 
      cand_e_size(1:nthreads+1) = 0
      address_cand_n(1:nthreads+1) = 0 
      address_cand_e(1:nthreads+1) = 0
 
      DO i=1,nsn
        iix(i)=0
        iiy(i)=0
        iiz(i)=0
        IF(stfn(i) == zero)cycle
        j=nsv(i)
 
C First test vs CRVOXEL
        ix=int(lrvoxel*(x(1,j)-xmin)/(xmax-xmin))
        IF(ix < 0 .OR. ix > lrvoxel) cycle
        iy=int(lrvoxel*(x(2,j)-ymin)/(ymax-ymin))
        IF(iy < 0 .OR. iy > lrvoxel) cycle
        iz=int(lrvoxel*(x(3,j)-zmin)/(zmax-zmin))
        IF(iz < 0 .OR. iz > lrvoxel) cycle
        IF(.NOT.(btest(crvoxel(iy,iz),ix))) cycle
 
C Second test vs reduced Box
        IF( (x(1,j)-xminb)/(xmaxb-xminb) > one ) THEN
          iix(i) = nbx
        ELSE
          iix(i)=int(max(nbx*(x(1,j)-xminb)/(xmaxb-xminb),-one))
        ENDIF
        IF( (x(2,j)-yminb)/(ymaxb-yminb) > one ) THEN
          iiy(i) = nby
        ELSE
          iiy(i)=int(max(nby*(x(2,j)-yminb)/(ymaxb-yminb),-one))
        ENDIF
        IF( (x(3,j)-zminb)/(zmaxb-zminb) > one ) THEN
          iiz(i) = nbz
        ELSE
          iiz(i)=int(max(nbz*(x(3,j)-zminb)/(zmaxb-zminb),-one))
        ENDIF
 
 
        iix(i)=max(1,2+iix(i))
        iiy(i)=max(1,2+iiy(i))
        iiz(i)=max(1,2+iiz(i))
 
        first = voxel(iix(i),iiy(i),iiz(i))
        IF(first == 0)THEN
c         Empty Cell
          voxel(iix(i),iiy(i),iiz(i)) = i ! first
          local_next_nod(i)                 = 0 ! last one
          last_nod(i)                 = 0 ! no last
        ELSEIF(last_nod(first) == 0)THEN
c         cell containing one node
c         add as next node
          local_next_nod(first) = i ! next
          last_nod(first) = i ! last
          local_next_nod(i)     = 0 ! last one
        ELSE
c
c         jump to the last node of the cell
          last = last_nod(first) ! last node in this voxel
          local_next_nod(last)  = i ! next
          last_nod(first) = i ! last
          local_next_nod(i)     = 0 ! last one
        ENDIF
      ENDDO
!$OMP END MASTER
 
!$OMP BARRIER
C=======================================================================
C 3 Searching boxes concerning each facet
C     and creation of candidates
C=======================================================================
 
      stagnod = numnod+numfakenodigeo
      IF(is_used_with_law151) stagnod = stagnod + numels
      ALLOCATE(tagnod(stagnod)) ; tagnod(1:stagnod) = 0
 
      local_i_stok = 0
      j_stok = 0
!$OMP DO SCHEDULE(guided)
      DO kk=1,nrtm_l
        ne=index(kk)
        IF(stf(ne) == zero)cycle
 
        IF(flagremnode==2.AND.iremnode==2)THEN
          k = kremnode(ne)+1
          l = kremnode(ne+1)
          DO m=k,l
            tagnod(remnode(m)) = 1
          ENDDO
        ENDIF
 
        IF(igap == 0)THEN
          aaa = tzinf+curv_max(ne)
        ELSE
          aaa = marge+curv_max(ne)+max(min(gapmax,max(gapmin,bgapsmx+gap_m(ne)))+dgapload,drad)
        ENDIF
 
        m1 = irect(1,ne)
        m2 = irect(2,ne)
        m3 = irect(3,ne)
        m4 = irect(4,ne)
 
        xx1=x(1,m1)
        xx2=x(1,m2)
        xx3=x(1,m3)
        xx4=x(1,m4)
        xmaxe=max(xx1,xx2,xx3,xx4)
        xmine=min(xx1,xx2,xx3,xx4)
 
        yy1=x(2,m1)
        yy2=x(2,m2)
        yy3=x(2,m3)
        yy4=x(2,m4)
        ymaxe=max(yy1,yy2,yy3,yy4)
        ymine=min(yy1,yy2,yy3,yy4)
 
        zz1=x(3,m1)
        zz2=x(3,m2)
        zz3=x(3,m3)
        zz4=x(3,m4)
        zmaxe=max(zz1,zz2,zz3,zz4)
        zmine=min(zz1,zz2,zz3,zz4)
 
 
c        calculation of the surface (for future elimination of candidates)
 
        sx = (yy3-yy1)*(zz4-zz2) - (zz3-zz1)*(yy4-yy2)
        sy = (zz3-zz1)*(xx4-xx2) - (xx3-xx1)*(zz4-zz2)
        sz = (xx3-xx1)*(yy4-yy2) - (yy3-yy1)*(xx4-xx2)
        s2 = sx*sx + sy*sy + sz*sz
 
c        index of voxels occupied by the facet
        IF( (xmine - aaa - xminb)/(xmaxb-xminb) > one ) THEN
          ix1 = nbx
        ELSE
          ix1=int(max(nbx*(xmine-aaa-xminb)/(xmaxb-xminb),-one))
        ENDIF
        IF( (ymine - aaa - yminb)/(ymaxb-yminb) > one ) THEN
          iy1 = nby
        ELSE
          iy1=int(max(nby*(ymine-aaa-yminb)/(ymaxb-yminb),-one))
        ENDIF
        IF( (zmine - aaa - zminb)/(zmaxb-zminb) > one ) THEN
          iz1 = nbz
        ELSE
          iz1=int(max(nbz*(zmine-aaa-zminb)/(zmaxb-zminb),-one))
        ENDIF
 
        ix1=max(1,2+ix1)
        iy1=max(1,2+iy1)
        iz1=max(1,2+iz1)
 
        IF( (xmaxe + aaa - xminb)/(xmaxb-xminb) > one ) THEN
          ix2 = nbx
        ELSE
          ix2=int(max(nbx*(xmaxe+aaa-xminb)/(xmaxb-xminb),-one))
        ENDIF
        IF( (ymaxe + aaa - yminb)/(ymaxb-yminb) > one ) THEN
          iy2 = nby
        ELSE
          iy2=int(max(nby*(ymaxe+aaa-yminb)/(ymaxb-yminb),-one))
        ENDIF
        IF( (zmaxe + aaa - zminb)/(zmaxb-zminb) > one ) THEN
          iz2 = nbz
        ELSE
          iz2=int(max(nbz*(zmaxe+aaa-zminb)/(zmaxb-zminb),-one))
        ENDIF
 
        ix2=max(1,2+ix2)
        iy2=max(1,2+iy2)
        iz2=max(1,2+iz2)
 
        DO iz = iz1,iz2
          DO iy = iy1,iy2
            DO ix = ix1,ix2
 
              jj = voxel(ix,iy,iz)
 
              DO WHILE(jj /= 0)
 
                nn=nsv(jj)
                IF(tagnod(nn) == 1 ) GOTO 300
                IF(nn == m1)GOTO 300
                IF(nn == m2)GOTO 300
                IF(nn == m3)GOTO 300
                IF(nn == m4)GOTO 300
 
                xs = x(1,nn)
                ys = x(2,nn)
                zs = x(3,nn)
                IF(igap /= 0)THEN
                  aaa = marge+curv_max(ne)+max(min(gapmax,max(gapmin,gap_s(jj)+gap_m(ne)))+dgapload,drad)
                ENDIF
                IF(xs<=xmine-aaa)GOTO 300
                IF(xs>=xmaxe+aaa)GOTO 300
                IF(ys<=ymine-aaa)GOTO 300
                IF(ys>=ymaxe+aaa)GOTO 300
                IF(zs<=zmine-aaa)GOTO 300
                IF(zs>=zmaxe+aaa)GOTO 300
 
c    underestimation of the distance**2 for candidate elimination
 
                d1x = xs - xx1
                d1y = ys - yy1
                d1z = zs - zz1
                d2x = xs - xx2
                d2y = ys - yy2
                d2z = zs - zz2
                dd1 = d1x*sx+d1y*sy+d1z*sz
                dd2 = d2x*sx+d2y*sy+d2z*sz
                IF(dd1*dd2 > zero)THEN
                  d2 = min(dd1*dd1,dd2*dd2)
                  a2 = aaa*aaa*s2
                  IF(d2 > a2)GOTO 300
                ENDIF
                j_stok = j_stok + 1
                prov_n(j_stok) = jj
                prov_e(j_stok) = ne
                IF(j_stok == nvsiz) THEN
                  CALL i7cor3(x    ,irect,nsv  ,prov_e ,prov_n,
     .                        stf  ,stfn ,gapv ,igap   ,gap   ,
     .                        gap_s,gap_m,istf ,gapmin ,gapmax,
     .                        gap_s_l,gap_m_l  ,drad   ,ix11,ix12,
     5                        ix13    ,ix14    ,nsvg,x1     ,x2     ,
     6                        x3     ,x4     ,y1  ,y2      ,y3      ,
     7                        y4     ,z1     ,z2  ,z3      ,z4      ,
     8                        xi     ,yi     ,zi  ,stif    ,dgapload,
     9                        j_stok)
                  CALL i7dst3(ix13,ix14,x1 ,x2 ,x3 ,
     1                 x4 ,y1 ,y2 ,y3 ,y4 ,
     2                 z1 ,z2 ,z3 ,z4 ,xi ,
     3                 yi ,zi ,x0 ,y0 ,z0 ,
     4                 nx1,ny1,nz1,nx2,ny2,
     5                 nz2,nx3,ny3,nz3,nx4,
     6                 ny4,nz4,p1 ,p2 ,p3 ,
     7                 p4 ,lb1,lb2,lb3,lb4,
     8                 lc1,lc2,lc3,lc4,j_stok)
                  CALL i7pen3(marge,gapv,n11,n21,n31,
     1                 pene ,nx1 ,ny1,nz1,nx2,
     2                 ny2  ,nz2 ,nx3,ny3,nz3,
     3                 nx4  ,ny4 ,nz4,p1 ,p2 ,
     4                 p3   ,p4,j_stok)
                  call inter_save_candidate( local_i_stok,j_stok,prov_n,prov_e,pene,local_cand_n,local_cand_e )
                  j_stok = 0
                ENDIF
 
  300           CONTINUE
                jj = local_next_nod(jj)
              ENDDO ! WHILE(JJ /= 0)
            ENDDO
          ENDDO
        ENDDO
 
        IF(flagremnode==2.AND.iremnode==2)THEN
          k = kremnode(ne)+1
          l = kremnode(ne+1)
          DO m=k,l
            tagnod(remnode(m)) = 0
          ENDDO
        ENDIF
 
      ENDDO
!$OMP END DO
      DEALLOCATE(tagnod)
 
      IF(j_stok/=0) THEN
        CALL i7cor3(x    ,irect,nsv  ,prov_e ,prov_n,
     .              stf  ,stfn ,gapv ,igap   ,gap   ,
     .              gap_s,gap_m,istf ,gapmin ,gapmax,
     .              gap_s_l,gap_m_l  ,drad   ,ix11,ix12,
     5              ix13    ,ix14    ,nsvg,x1     ,x2     ,
     6              x3     ,x4     ,y1  ,y2      ,y3      ,
     7              y4     ,z1     ,z2  ,z3      ,z4      ,
     8              xi     ,yi     ,zi  ,stif    ,dgapload,
     9              j_stok)
        CALL i7dst3(ix13,ix14,x1 ,x2 ,x3 ,
     1              x4 ,y1 ,y2 ,y3 ,y4 ,
     2              z1 ,z2 ,z3 ,z4 ,xi ,
     3              yi ,zi ,x0 ,y0 ,z0 ,
     4              nx1,ny1,nz1,nx2,ny2,
     5              nz2,nx3,ny3,nz3,nx4,
     6              ny4,nz4,p1 ,p2 ,p3 ,
     7              p4 ,lb1,lb2,lb3,lb4,
     8              lc1,lc2,lc3,lc4,j_stok)
        CALL i7pen3(marge,gapv,n11,n21,n31,
     1              pene ,nx1 ,ny1,nz1,nx2,
     2              ny2  ,nz2 ,nx3,ny3,nz3,
     3              nx4  ,ny4 ,nz4,p1 ,p2 ,
     4              p3   ,p4,j_stok)
        call inter_save_candidate( local_i_stok,j_stok,prov_n,prov_e,pene,local_cand_n,local_cand_e )
        j_stok = 0
      ENDIF
 
      ! save the local number of candidates
      cand_n_size(itask+1) = local_i_stok
      cand_e_size(itask+1) = local_i_stok
!$OMP BARRIER
 
 
!$OMP SINGLE
      ! compute the address for each threads & the total number of candidates
      address_cand_n(1) = 0
      address_cand_e(1) = 0
      ! ------------
      do i=1,nthreads
        address_cand_n(i+1) = cand_n_size(i) + address_cand_n(i) ! address for cand_n
        address_cand_e(i+1) = cand_e_size(i) + address_cand_e(i) ! address for cand_e
 
        cand_n_size(nthreads+1) = cand_n_size(nthreads+1) + cand_n_size(i) ! total number of cand_n for this proc
        cand_e_size(nthreads+1) = cand_e_size(nthreads+1) + cand_e_size(i) ! total number of cand_e for this proc
      enddo
      ! ------------
 
      ! ------------
      ! check the size of global cand_n & cand_e
      my_old_size = ipari(18)*ipari(23)
      i_stok_old = i_stok 
      i_stok =  i_stok + cand_n_size(nthreads+1) ! total number of cand_n/cand_e (= sum( cand_n/e per proc))
      if(i_stok > my_old_size) then ! total number of candidates > size of cand_n/e --> need to extend the 2 arrays
        multimp = i_stok/ipari(18) + 1
        call upgrade_multimp(nin,multimp,intbuf_tab)
      endif
      ! ------------
!$OMP END SINGLE
 
!$OMP BARRIER
 
      ! ------------
      my_address = i_stok_old + address_cand_n(itask+1) ! get the address for each thread
      intbuf_tab%cand_n(my_address+1:my_address+local_i_stok) = local_cand_n%int_array_1d(1:local_i_stok) ! save the cand_n into the global array
      my_address = i_stok_old + address_cand_e(itask+1) ! get the address for each thread
      intbuf_tab%cand_e(my_address+1:my_address+local_i_stok) = local_cand_e%int_array_1d(1:local_i_stok) ! save the cand_e into the global array
      ! ------------
 
 
      call dealloc_1d_array(local_cand_n)
      call dealloc_1d_array(local_cand_e)
 
!$OMP BARRIER
 
      ! Sort the candidates to ensure the same domain decomposition
!$OMP SINGLE
      ! ------------
      my_size = cand_n_size(nthreads+1)
      allocate(my_index(2*my_size))
      allocate(sort_array(2,my_size))
      allocate(save_array(2,my_size))
      
      my_address = i_stok_old + address_cand_n(1) ! get the address of the first pair of candidate
      sort_array(1,1:my_size) = intbuf_tab%cand_n(my_address+1:my_address+my_size)
      my_address = i_stok_old + address_cand_e(1) ! get the address of the first pair of candidate
      sort_array(2,1:my_size) = intbuf_tab%cand_e(my_address+1:my_address+my_size)
      do i=1,my_size
        my_index(i) = i
      enddo
      save_array(1:2,1:my_size) = sort_array(1:2,1:my_size)
      mode = 0
 
      call my_orders( mode,work,sort_array,my_index,my_size,2)
      my_address = i_stok_old + address_cand_n(1) ! get the address of the first pair of candidate
      do i=1,my_size 
        intbuf_tab%cand_n(my_address+i) = save_array(1,my_index(i))
      enddo 
      my_address = i_stok_old + address_cand_e(1) ! get the address of the first pair of candidate
      do i=1,my_size
        intbuf_tab%cand_e(my_address+i) = save_array(2,my_index(i))
      enddo
      deallocate(my_index)
      deallocate(sort_array)
      deallocate(save_array)
      ! ------------
!$OMP END SINGLE    
 
!$OMP DO SCHEDULE(guided)
      DO i=1,nsn
        IF(iix(i)/=0)THEN
          voxel(iix(i),iiy(i),iiz(i))=0
        ENDIF
      ENDDO
!$OMP END DO
 
!$OMP SINGLE
      deallocate( cand_n_size,cand_e_size )
      deallocate( address_cand_n,address_cand_e )
!$OMP END SINGLE
 
      RETURN
      END