 SUBROUTINE block2mesh(cb,nzone,zone,nzonei,zonei)
!=========================== DOCUMENTATION ============================= 
!DOCDEBh*    UTI/PRE/block2mesh  *
!
!    NOM
! block2mesh(cb,nzone,zone,nzonei,zonei)
!    DESCRIPTION  
!*
! A partir du maillage maitre (les blocs) on construit le maillage de chaque domaine.
! From the master mesh (the blocks) we construct the mesh for each domain.
!*
!
! 
! 
!    SOURCE
! Ce logiciel est regi par la licence [CeCILL-B]
! This software is governed by the [CeCILL-B] license
!=========================== DEBUT DES DECLARATIONS ==================== 
!.1-----  Implicit, Use
  use num
  use uti, only: inter,init_def_bloc,ini_w,splitdom_Z ! ,bathy_sv1,initsu_sv1
  IMPLICIT NONE
!.2-----  Declaration
  type(calcul)                                      :: cb !! Objet calcul contenant le maillage maitre / Calcul object containing the master mesh
  integer                                           :: nzonei  
!! Parametres de definition des maillages et conditions aux limites / Mesh definition parameters and boundary conditions
  real(kind=kind(0.d0)),allocatable,dimension(:,:)  :: zone,zonei
!! Parametres de definition des maillages et conditions aux limites / Mesh definition parameters and boundary conditions
  integer,allocatable,dimension(:,:)                :: nzone
!! Parametres de definition des maillages et conditions aux limites / Mesh definition parameters and boundary conditions

  type(bloc)                                        :: blv(6),bl_null
  type(calcul)                                      :: calc
  integer                                           :: imac(6),imav(6),imaf(6)
  integer                                           :: idom,nb_dom,nx,ny,nz,i,nb_bloc,ib
  integer                                           :: icb,ibv
  integer                                           :: k,nvi
  integer                                           :: l
  real(kind=kind(0.d0))                             :: hmin,dd,x,y
  type(list_int), allocatable,dimension(:)          :: interdom
  integer, allocatable,dimension(:)                 :: nbl

!=========================== DEBUT DU CODE EXECUTABLE ==================
bl_null%numbl =0 ;bl_null%idom =0 ;bl_null%morton  =0 ;bl_null%nrb =0 ;bl_null%nx =0 ;bl_null%ny =0 ;bl_null%nz  =0 ;
bl_null%volc  =0.d0 ;bl_null%cl(:)  =0 ;bl_null%nrf(:)  =0 ;bl_null%dom(:)  =0 ;bl_null%first_cell =0 ;
bl_null%nb_cell   =0 ;   bl_null%first_face =0 ;bl_null%nb_face    =0 ;  bl_null%first_vertex =0 ;bl_null%nb_vertex  =0 


  write(*,*)
  write(*,*)achar(27)//'[34m================================================'
  write(*,*)'Creation of the domain mesh from the master mesh'
  write(*,*)'================================================'//achar(27)//'[0m'
  write(*,*)
  nb_bloc=cb%nb_bloc
  nb_dom=cb%nbds
  write(*,*)'BLOCK2MESH-INFO Problem dimension :',cb%msh%ndim
! estimation de hmin
    hmin=1.d20
    DO ib=1,nb_bloc
        dd=sqrt((cb%list_bloc(ib)%s(1)%x-cb%list_bloc(ib)%s(2)%x)**2+ &
                (cb%list_bloc(ib)%s(1)%y-cb%list_bloc(ib)%s(2)%y)**2+(cb%list_bloc(ib)%s(1)%z-cb%list_bloc(ib)%s(2)%z)**2)
        if(dd.lt.hmin)hmin=dd
        dd=sqrt((cb%list_bloc(ib)%s(1)%x-cb%list_bloc(ib)%s(4)%x)**2+ &
                (cb%list_bloc(ib)%s(1)%y-cb%list_bloc(ib)%s(4)%y)**2+(cb%list_bloc(ib)%s(1)%z-cb%list_bloc(ib)%s(4)%z)**2)
        if(dd.lt.hmin)hmin=dd
        dd=sqrt((cb%list_bloc(ib)%s(1)%x-cb%list_bloc(ib)%s(5)%x)**2+ &
                (cb%list_bloc(ib)%s(1)%y-cb%list_bloc(ib)%s(5)%y)**2+(cb%list_bloc(ib)%s(1)%z-cb%list_bloc(ib)%s(5)%z)**2)
        if(dd.lt.hmin)hmin=dd
    END DO

!----------------------------------------------------------------------------------
! on definit la discretisation des blocs et l appartenance a un domaine
!
  CALL init_def_bloc(cb,nzone,zone)

!----------------------------------------------------------------------------------
! on verifie la condition sur le niveau de discretisation de 2 blocs voisins
  DO l=1,cb%nrma
    DO ib=1,nb_bloc
        DO i=1,6
            IF (cb%list_bloc(ib)%cl(i).gt.0) THEN
               ibv=cb%list_bloc(ib)%cl(i)
               IF(cb%list_bloc(ib)%nrb.eq.(cb%nrma+1-l).and.cb%list_bloc(ibv)%nrb.lt.(cb%nrma-l)) cb%list_bloc(ibv)%nrb=(cb%nrma-l)
            END IF
        END DO
    END DO
  END DO
!----------------------------------------------------------------------------------
! on decompose en sous-domaines par code de Morton

  write( *,1011)nb_bloc,nb_dom
1011 format('Aggregation of ',i6,' blocks into',i3,' sub-domains')

  CALL splitdom_Z(cb%list_bloc,nb_bloc,nb_dom, cb%msh%ndim) 

! liste des voisins: -1/ le voisin est + grossier
!                     0/ le voisin est identique
!                     1/ le voisin est + fin
  DO ib=1,nb_bloc
      cb%list_bloc(ib)%nrf=0
      DO i=1,6
          ibv=cb%list_bloc(ib)%cl(i)
          IF(ibv.lt.-100000)ibv=-ibv-100000
          IF (ibv.gt.0) THEN
              IF (cb%list_bloc(ibv)%nrb.gt.cb%list_bloc(ib)%nrb) cb%list_bloc(ib)%nrf(i)=1
              IF (cb%list_bloc(ibv)%nrb.lt.cb%list_bloc(ib)%nrb) cb%list_bloc(ib)%nrf(i)=-1
          END IF
      END DO
  END DO
!
!-------  Determination de la taille des tableaux du bloc
  DO ib=1,nb_bloc
      cb%list_bloc(ib)%dom=0
      DO i=1,6
          ibv=cb%list_bloc(ib)%cl(i)
          IF(ibv.lt.-100000)ibv=-ibv-100000
          IF (ibv.gt.0) cb%list_bloc(ib)%dom(i)=cb%list_bloc(ibv)%idom
      END DO
      imac=1
      imav=1
      imaf=1
      DO i=1,6 ! ima: indice de maillage de l interface entre 2 blocs voisins
          ibv=cb%list_bloc(ib)%cl(i)
          IF(ibv.lt.-100000)ibv=-ibv-100000
          IF (cb%list_bloc(ib)%idom.eq.cb%list_bloc(ib)%dom(i).and.ib.lt.ibv) imaf(i)=0
          IF (cb%list_bloc(ib)%idom.eq.cb%list_bloc(ib)%dom(i).and.ib.lt.ibv) imav(i)=0
          IF (cb%list_bloc(ib)%idom.eq.cb%list_bloc(ib)%dom(i).or.ibv.le.0)   imac(i)=0
      END DO

      SELECT CASE (cb%msh%ndim)
        CASE (1) ! Cas 1D
          nx=cb%list_bloc(ib)%nx*2**(cb%list_bloc(ib)%nrb)
          ny=1
          nz=1
          cb%list_bloc(ib)%nb_cell   = nx+imac(1)+imac(2)
          cb%list_bloc(ib)%nb_vertex = (nx+1)*(ny+1)*2 
          cb%list_bloc(ib)%nb_face   = (nx-1)+nx*2*2 + imaf(1)+imaf(2) 

         CASE (2) ! Cas 2D
          nx=cb%list_bloc(ib)%nx*2**(cb%list_bloc(ib)%nrb)
          ny=cb%list_bloc(ib)%ny*2**(cb%list_bloc(ib)%nrb)
          nz=1

          cb%list_bloc(ib)%nb_cell   = nx*ny + int(ny*2.**cb%list_bloc(ib)%nrf(1)*imac(1)) &
                                   + int(ny*2.**cb%list_bloc(ib)%nrf(2)*imac(2)) &
                                   + int(nx*2.**cb%list_bloc(ib)%nrf(3)*imac(3)) &
                                   + int(nx*2.**cb%list_bloc(ib)%nrf(4)*imac(4))
          cb%list_bloc(ib)%nb_vertex = (nx+1)*(ny+1)*2 &
                                   + cb%list_bloc(ib)%nrf(1)*(cb%list_bloc(ib)%nrf(1)+1)*ny*imav(1) &
                                   + cb%list_bloc(ib)%nrf(2)*(cb%list_bloc(ib)%nrf(2)+1)*ny*imav(2) &
                                   + cb%list_bloc(ib)%nrf(3)*(cb%list_bloc(ib)%nrf(3)+1)*nx*imav(3) &
                                   + cb%list_bloc(ib)%nrf(4)*(cb%list_bloc(ib)%nrf(4)+1)*nx*imav(4)
          cb%list_bloc(ib)%nb_face   = (nx-1)*ny+nx*(ny-1)+nx*ny*2 &
                                   + ny*imaf(1)+ny*imaf(2)+nx*imaf(3)+nx*imaf(4)  &
                                   + cb%list_bloc(ib)%nrf(1)*(cb%list_bloc(ib)%nrf(1)+1)*ny/2*imaf(1) &
                                   + cb%list_bloc(ib)%nrf(2)*(cb%list_bloc(ib)%nrf(2)+1)*ny/2*imaf(2) &
                                   + cb%list_bloc(ib)%nrf(3)*(cb%list_bloc(ib)%nrf(3)+1)*nx/2*imaf(3) &
                                   + cb%list_bloc(ib)%nrf(4)*(cb%list_bloc(ib)%nrf(4)+1)*nx/2*imaf(4)        
         CASE (3) ! Cas 3D
          nx=cb%list_bloc(ib)%nx*2**(cb%list_bloc(ib)%nrb)
          ny=cb%list_bloc(ib)%ny*2**(cb%list_bloc(ib)%nrb)
          nz=cb%list_bloc(ib)%nz*2**(cb%list_bloc(ib)%nrb)
          cb%list_bloc(ib)%nb_cell   = nx*ny*nz &
                                   + int(ny*nz*4.**cb%list_bloc(ib)%nrf(1)*imac(1)) &
                                   + int(ny*nz*4.**cb%list_bloc(ib)%nrf(2)*imac(2)) &
                                   + int(nx*nz*4.**cb%list_bloc(ib)%nrf(3)*imac(3)) &
                                   + int(nx*nz*4.**cb%list_bloc(ib)%nrf(4)*imac(4)) &
                                   + int(nx*ny*4.**cb%list_bloc(ib)%nrf(5)*imac(5)) &
                                   + int(nx*ny*4.**cb%list_bloc(ib)%nrf(6)*imac(6)) 
          cb%list_bloc(ib)%nb_vertex = (nx+1)*(ny+1)*(nz+1) &
                                   + cb%list_bloc(ib)%nrf(1)*(cb%list_bloc(ib)%nrf(1)+1)/2*(3*ny*nz+ny+nz)*imav(1) &
                                   + cb%list_bloc(ib)%nrf(2)*(cb%list_bloc(ib)%nrf(2)+1)/2*(3*ny*nz+ny+nz)*imav(2) &
                                   + cb%list_bloc(ib)%nrf(3)*(cb%list_bloc(ib)%nrf(3)+1)/2*(3*nx*nz+nx+nz)*imav(3) &
                                   + cb%list_bloc(ib)%nrf(4)*(cb%list_bloc(ib)%nrf(4)+1)/2*(3*nx*nz+nx+nz)*imav(4) &
                                   + cb%list_bloc(ib)%nrf(5)*(cb%list_bloc(ib)%nrf(5)+1)/2*(3*nx*ny+nx+ny)*imav(5) &
                                   + cb%list_bloc(ib)%nrf(6)*(cb%list_bloc(ib)%nrf(6)+1)/2*(3*nx*ny+nx+ny)*imav(6) 
          cb%list_bloc(ib)%nb_face   = (nx-1)*ny*nz+nx*(ny-1)*nz+nx*ny*(nz-1) &
                                   + ny*nz*imaf(1)+ny*nz*imaf(2)+nx*nz*imaf(3)+nx*nz*imaf(4) &
                                   + nx*ny*imaf(5)+nx*ny*imaf(6) &
                                   + cb%list_bloc(ib)%nrf(1)*(cb%list_bloc(ib)%nrf(1)+1)/2*3*ny*nz*imaf(1) &
                                   + cb%list_bloc(ib)%nrf(2)*(cb%list_bloc(ib)%nrf(2)+1)/2*3*ny*nz*imaf(2) &
                                   + cb%list_bloc(ib)%nrf(3)*(cb%list_bloc(ib)%nrf(3)+1)/2*3*nx*nz*imaf(3) &
                                   + cb%list_bloc(ib)%nrf(4)*(cb%list_bloc(ib)%nrf(4)+1)/2*3*nx*nz*imaf(4) &  
                                   + cb%list_bloc(ib)%nrf(5)*(cb%list_bloc(ib)%nrf(5)+1)/2*3*nx*ny*imaf(5) &  
                                   + cb%list_bloc(ib)%nrf(6)*(cb%list_bloc(ib)%nrf(6)+1)/2*3*nx*ny*imaf(6)   
                                                      
      END SELECT
      !call display(cb%list_bloc(ib)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  END DO

! Creation de la liste des interfaces 
  allocate(interdom(int(nb_dom*(nb_dom-1)/2)))
  allocate(nbl(int(nb_dom*(nb_dom-1)/2)))
  CALL inter(cb%list_bloc,nb_bloc,nb_dom,interdom,nbl,cb%msh%ndim)
  
! Determination de petit
  petit=1.d-3
  DO ib=1,nb_bloc
      IF (hmin/(cb%list_bloc(ib)%nx*2**cb%list_bloc(ib)%nrb).lt.petit) &
          petit=hmin/(cb%list_bloc(ib)%nx*2**cb%list_bloc(ib)%nrb)
      IF (hmin/(cb%list_bloc(ib)%ny*2**cb%list_bloc(ib)%nrb).lt.petit.and.cb%msh%ndim.ne.1) &
          petit=hmin/(cb%list_bloc(ib)%ny*2**cb%list_bloc(ib)%nrb)
      IF (hmin/(cb%list_bloc(ib)%nz*2**cb%list_bloc(ib)%nrb).lt.petit.and.cb%msh%ndim.eq.3) &
          petit=hmin/(cb%list_bloc(ib)%nz*2**cb%list_bloc(ib)%nrb)
  END DO
  petit=petit*1.d-5
  write( *,1212)petit
 1212  format('Small dimension',e13.5)

 !----------- Cas particulier de Saint Venant -----
! On stocke eventuellement la bathy sur chaque bloc
if(imodel_phy.eq.1.and.ibathy.eq.2)then
open(14,file=fich_bathy,status='old')
   do ib=1,nb_bloc
       do i=1,4
          read(14,*)x,y,cb%list_bloc(ib)%vb(i), cb%list_bloc(ib)%vb(i+4)
       enddo
   enddo
   close(14)
endif


!-----------------------------------------------------------------------
! Boucle sur les domaines pour creer calc
!-----------------------------------------------------------------------
  DO idom=1,nb_dom
! Initialisations
      calc%msh%ndim  =cb%msh%ndim
      calc%msh%nb_dom=nb_dom
      nullify(calc%msh%list_vertex)
      nullify(calc%msh%list_face)
      nullify(calc%msh%list_cell)
      nullify(calc%msh%list_send)
      nullify(calc%msh%list_recv)
      CALL empty_mesh(calc%msh,calc%msh%nb_dom,idom)
      IF (associated(calc%list_bloc)) THEN
          deallocate(calc%list_bloc)
          nullify(calc%list_bloc)
      END IF
      nullify(calc%list_bloc)
     
! Determination de la taille des tableaux du domaine
      calc%msh%nb_cell   = 0
      calc%msh%nb_vertex = 0
      calc%msh%nb_face   = 0
      calc%nb_bloc       = 0
         
      DO ib=1,nb_bloc
          IF (cb%list_bloc(ib)%idom.eq.idom) THEN
              calc%nb_bloc=calc%nb_bloc+1
              cb%list_bloc(ib)%first_cell   = calc%msh%nb_cell
              calc%msh%nb_cell    = calc%msh%nb_cell + cb%list_bloc(ib)%nb_cell
              cb%list_bloc(ib)%first_face   = calc%msh%nb_face
              calc%msh%nb_face    = calc%msh%nb_face + cb%list_bloc(ib)%nb_face
              cb%list_bloc(ib)%first_vertex = calc%msh%nb_vertex
              calc%msh%nb_vertex  = calc%msh%nb_vertex + cb%list_bloc(ib)%nb_vertex
          END IF
      END DO
  
      allocate(calc%list_bloc(calc%nb_bloc))
      allocate(calc%msh%list_cell(calc%msh%nb_cell))
      allocate(calc%msh%list_vertex(calc%msh%nb_vertex))
      allocate(calc%msh%list_face(calc%msh%nb_face))
      allocate(calc%msh%list_send(calc%msh%nb_dom))
      allocate(calc%msh%list_recv(calc%msh%nb_dom))
  
      DO i=1,idom-1
          k=idom-i+(i-1)*(2*nb_dom-i)/2 
          IF (nbl(k).ne.0) allocate(calc%msh%list_recv(i)%L(interdom(k)%L(4*nbl(k)-1)))
          IF (nbl(k).ne.0) allocate(calc%msh%list_send(i)%L(interdom(k)%L(4*nbl(k)-0)))
      END DO
      nullify(calc%msh%list_send(idom)%L)
      nullify(calc%msh%list_recv(idom)%L)
      DO i=idom+1,nb_dom
          k=i-idom+(idom-1)*(2*nb_dom-idom)/2
          IF (nbl(k).ne.0) allocate(calc%msh%list_recv(i)%L(interdom(k)%L(4*nbl(k)-0)))
          IF (nbl(k).ne.0) allocate(calc%msh%list_send(i)%L(interdom(k)%L(4*nbl(k)-1)))
      END DO
      CALL initialise(calc%msh)
       
! Maillage des blocs du domaine
      icb=0
      DO ib=1,nb_bloc
          IF (cb%list_bloc(ib)%idom.eq.idom) THEN
              icb=icb+1
              calc%list_bloc(icb)=cb%list_bloc(ib)
              DO i=1,6
                  blv(i)=bl_null
                  ibv=cb%list_bloc(ib)%cl(i)
                  if(ibv.lt.-100000)ibv=-ibv-100000
                  IF (ibv.gt.0) blv(i)=cb%list_bloc(ibv)
              END DO

              SELECT CASE (cb%msh%ndim)
                 CASE (1)
                   CALL mesh1d(cb%list_bloc(ib),calc%msh,ib,nb_dom,blv,interdom)

                 CASE (2)
                   CALL mesh2d(cb%list_bloc(ib),calc%msh,ib,nb_dom,blv,interdom)

                 CASE (3)
                   CALL mesh3d(cb%list_bloc(ib),calc%msh,ib,nb_dom,blv,interdom)

            END SELECT
          END IF
      END DO ! Fin boucle sur les blocs

      CALL init_mesh(calc%msh)

! On elimine les sommets surabondants
      nvi=calc%msh%nb_vertex
      CALL trivertex(calc%msh)
      write( *,1111) idom,nvi-calc%msh%nb_vertex,nvi,calc%msh%nb_vertex
1111  format('Domain',i3,' elimination of ',i7,' vertexes out of',i7,', i.e.',i7,' vertexes.')
      write( *,2222) idom,calc%msh%nb_vertex
      write( *,3333) idom,calc%msh%nb_face
      write( *,4444) idom,calc%msh%nb_cell
2222  format('Domain',i3,' Number of vertexes :',i7)
3333  format('Domain',i3,' Number of faces    :',i7)
4444  format('Domain',i3,' Number of cells    :',i7)
        
! Recuperation des données de cb
      calc%tmin           = cb%tmin
      calc%tmax           = cb%tmax
      calc%dt             = cb%dt
      calc%pdtp           = cb%pdtp
      calc%cfl            = cb%cfl
      calc%vmax           = cb%vmax
      calc%isauv          = cb%isauv
      calc%istart         = cb%istart
      calc%ipas           = cb%ipas
      calc%iordre_e       = cb%iordre_e
      calc%iordre_t       = cb%iordre_t
      calc%imeth_int_temps= cb%imeth_int_temps
      calc%point_sondes   = cb%point_sondes
      calc%nb_sondes      = cb%nb_sondes
      calc%num_point_min  = cb%num_point_min
      calc%ityp_flux      = cb%ityp_flux
      calc%ilimiteur      = cb%ilimiteur
      calc%nbds           = cb%nbds
      calc%nrma           = cb%nrma
      calc%vcde           = cb%vcde
      calc%vcra           = cb%vcra
      if(cb%nb_sol.ne.0)call copy_sol(cb,calc)
      
! Initialisation sur tout le domaine
  
      calc%msh%numdom=idom
      call ini_w(calc,nzonei,zonei)
       !call display(calc%msh)  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      CALL save_calc(calc)
      CALL empty_mesh(calc%msh,0,0)  ! liberation de la memoire

  END DO ! Fin boucle sur les domaines

  !deallocate(bl)

!===========================   FIN DE LA ROUTINE    ====================

END SUBROUTINE block2mesh
!DOCFIN 
