  PROGRAM cerf_amr
!=========================== DOCUMENTATION ============================= 
!DOCDEBh*    UTI/cerf_amr  *
!
!    NOM
! convert_mai
!
!    DESCRIPTION  
!*
! Programme permettant la lecture de fichiers binaires au format CERF
! puis de calculer un critere afin de raffiner ou deraffiner eventuellement
!
! Program for reading binary files in CERF format
! then calculates a criterion for refining or coarsening the mesh if necessary
!*
! 
!    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 uti
  IMPLICIT NONE
!.2-----  Declaration
  type(bloc),allocatable,dimension(:)             :: bl_ini,bl_out
  type(bloc)                                      :: blv(6),bl_null
  type(list_vec),allocatable,dimension(:)         :: list_w_bloc
  type(calcul)                                    :: calc_ini,calc_out,cb
  integer                                         :: nbdom,nrafmax,i,iv,idv,j,k,nvi,bcl,bclj,ios
  integer                                         :: itest, idom, nbtotcelli,nbtotcinti,nbtotcello,nbtotcinto,nb_bloc
  integer                                         :: nb_cell_bloc,ib,ibloc,ic1,icb,itemp,nx,ny,nz
  real(kind=kind(0.d0))                           :: crit_raf,crit_deraf,critmoy,voltot
  real(kind=kind(0.d0))                           :: delta,d
  real(kind=kind(0.d0)), dimension(:),allocatable :: critere,dtb
  integer, dimension(:),allocatable               :: im
  integer                                         :: imac(6),imav(6),imaf(6)
  type(list_int), allocatable,dimension(:)        :: interdom
  integer, allocatable,dimension(:)               :: nbl,isb
  character(len=256)                              :: iomsg
!=========================== DEBUT DU CODE EXECUTABLE ==================

  bl_null%idom=0
! lecture de l objet calcul d ou on va extraire des blocs
  CALL empty_mesh(cb%msh,1,0)
  CALL get_calc(cb)
  nrafmax   =cb%nrma ! Niveau maximal de raffinement
  crit_raf  =cb%vcra ! critere a partir duquel on raffine
  crit_deraf=cb%vcde ! critere en dessous duquel on deraffine
  write( *,1051)nrafmax
  write( *,1052)crit_raf
  write( *,1053)crit_deraf
 1051 format('CERF_AMR-Info  Maximal refinement level :',i5)
 1052 format('CERF_AMR-Info Mesh refinement parameter :',e13.5)
 1053 format('CERF_AMR-Info Mesh coarsening parameter :',e13.5)

! On normalise les coordonnees des sommet des blocs a 1/1000 de hpetit
  hpetit=1.d20
  DO i=1,cb%msh%nb_cell
      IF (cb%msh%list_cell(i)%h.lt.hpetit)hpetit=cb%msh%list_cell(i)%h
  END DO
  OPEN(11,file='cerfbin_000',status='old',form='unformatted',iostat=ios,iomsg=iomsg)
  if (ios /= 0) call print_err('CERF_AMR','Cannot open cerfbin_000: '//trim(iomsg))
  read(11,iostat=ios,iomsg=iomsg) nb_bloc,nbdom
  if (ios /= 0) call print_err('CERF_AMR','Cannot read header from cerfbin_000: '//trim(iomsg))


  CLOSE(11)
  allocate(list_w_bloc(nb_bloc))
  allocate(bl_ini(nb_bloc))
  allocate(bl_out(nb_bloc))
  allocate(critere(nb_bloc))
  allocate(dtb(nb_bloc))
  allocate(im(nb_bloc))
  allocate(isb(nb_bloc)) ! Indice de solide par bloc
  isb=0

  im=0
  dtb=1.d20
  critere=0.d0

!-----------------------------------------------------------------------------------
! Boucle sur les domaines et les blocs pour determiner le critere de raffinement
!-----------------------------------------------------------------------------------
  critmoy=0.d0   ! valeur moyenne du critere sur tous les domaines
  voltot=0.d0    ! voulume total 
  nbtotcelli=0   !n ombre total (sur tous les domaines) de cellules interieures et aux interfaces
  nbtotcinti=0   ! nombre total (sur tous les domaines) de cellules interieures 
  hpetit=1.d20     ! dimension du plus petit bloc
  DO idom=1,nbdom
      CALL empty_mesh(calc_ini%msh,nbdom,idom)
      CALL empty_sol(calc_ini%list_sol)
      calc_ini%msh%nb_dom=nbdom
      calc_ini%ipas=cb%ipas
! Recuperation du calcul et sauvegarde des variables primitives internes au bloc
      CALL get_calc(calc_ini)
 !!!!!!!!!!!!!
 ! on met a jour les objets "solide" dans cb a partir du domaine 1
  if(idom.eq.1.and.cb%nb_sol.ne.0)then
    do i=1,cb%nb_sol
        cb%list_sol(i)%imv    =calc_ini%list_sol(i)%imv
        cb%list_sol(i)%rho    =calc_ini%list_sol(i)%rho
        cb%list_sol(i)%xmin_s =calc_ini%list_sol(i)%xmin_s
        cb%list_sol(i)%xmax_s =calc_ini%list_sol(i)%xmax_s
        cb%list_sol(i)%ymin_s =calc_ini%list_sol(i)%ymin_s
        cb%list_sol(i)%ymax_s =calc_ini%list_sol(i)%ymax_s
        cb%list_sol(i)%zmin_s =calc_ini%list_sol(i)%zmin_s
        cb%list_sol(i)%zmax_s =calc_ini%list_sol(i)%zmax_s
        cb%list_sol(i)%xg     =calc_ini%list_sol(i)%xg
        cb%list_sol(i)%theta  =calc_ini%list_sol(i)%theta
        cb%list_sol(i)%list_vertex  =calc_ini%list_sol(i)%list_vertex
    enddo
  endif
 !!!!!!!!!!!!!
! Boucle sur les blocs du domaine
      DO ib=1,calc_ini%nb_bloc
          ibloc=calc_ini%list_bloc(ib)%numbl
          bl_ini(ibloc)=calc_ini%list_bloc(ib)
! determination du nombre de cellules du bloc et hpetit, dimension du plus petit bloc
          SELECT CASE (calc_ini%msh%ndim)
             CASE (1) ! Cas 1D
              nb_cell_bloc=bl_ini(ibloc)%nx*2**bl_ini(ibloc)%nrb
              IF (abs(bl_ini(ibloc)%s(1)%x-bl_ini(ibloc)%s(2)%x).lt.hpetit)&
              hpetit=abs(bl_ini(ibloc)%s(1)%x-bl_ini(ibloc)%s(2)%x)
             CASE (2) ! Cas 2D
              nb_cell_bloc=bl_ini(ibloc)%nx*2**bl_ini(ibloc)%nrb* &
                           bl_ini(ibloc)%ny*2**bl_ini(ibloc)%nrb
              IF (abs(bl_ini(ibloc)%s(1)%y-bl_ini(ibloc)%s(4)%y).lt.hpetit)&
              hpetit=abs(bl_ini(ibloc)%s(1)%y-bl_ini(ibloc)%s(4)%y)
             CASE (3) ! Cas 3D
              nb_cell_bloc=bl_ini(ibloc)%nx*2**bl_ini(ibloc)%nrb* &
                           bl_ini(ibloc)%ny*2**bl_ini(ibloc)%nrb* &
                           bl_ini(ibloc)%nz*2**bl_ini(ibloc)%nrb
              IF (abs(bl_ini(ibloc)%s(1)%z-bl_ini(ibloc)%s(5)%z).lt.hpetit)&
              hpetit=abs(bl_ini(ibloc)%s(1)%z-bl_ini(ibloc)%s(5)%z)
          END SELECT
! Sauvegarde des variables primitives internes au bloc
          nullify(list_w_bloc(ibloc)%vec)
          allocate(list_w_bloc(ibloc)%vec(nb_cell_bloc))
          DO i=1,nb_cell_bloc
              list_w_bloc(ibloc)%vec(i)=calc_ini%msh%list_cell(bl_ini(ibloc)%first_cell+i)%w%vprim
              if(calc_ini%msh%list_cell(bl_ini(ibloc)%first_cell+i)%w%chi.ne.0.d0)isb(ibloc)=1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          END DO
! Calcul du critere
          CALL critraf(calc_ini,critere,ib,nb_cell_bloc)
          critmoy=critmoy+critere(ibloc)
          voltot=voltot+volbloc(bl_ini(ibloc))
          nbtotcinti=nbtotcinti+nb_cell_bloc  
      END DO
      nbtotcelli=nbtotcelli+calc_ini%msh%nb_cell
      
! Reecriture sur fichier solbin
      delta = calc_ini%tmax-calc_ini%tmin
      calc_ini%tmin  =calc_ini%tmin-delta
      calc_ini%tmax  =calc_ini%tmax-delta           
      CALL save_calc(calc_ini)
      CALL empty_mesh(calc_ini%msh,0,0)  ! liberation de la mémoire
      CALL empty_sol(calc_ini%list_sol)
  END DO
!------------------------------------------------------------------------------------------
! FIN de la boucle sur les domaines et les blocs pour determiner le critere de raffinement
!-------------------------------------------------------------------------------------------
  do ib=1,nb_bloc
    critere(ib)=critere(ib)/volbloc(bl_ini(ib))
  enddo
  if(voltot.eq.0.d0)then
      call print_err('CERF_AMR','Total volume is zero, cannot compute average criterion')
  else
      critmoy=critmoy/voltot
  endif
  if(critmoy.eq.0.d0)then
      call print_err('CERF_AMR','Average criterion is zero, cannot normalize criterion')
  else
      critere=critere/critmoy
      write(*,*)'CERF_AMR-INFO Average criterion : ',critmoy
  endif

! Normalisation du critere
  CALL amr2tec(bl_ini,critere,nb_bloc) 
  bl_out=bl_ini
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  DO ib=1,nb_bloc
    if(isb(ib).eq.1)critere(ib)=2*crit_raf
  enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Determination des niveaux de raffinement des blocs
  DO ib=1,nb_bloc
      IF (dabs(critere(ib)).gt.crit_raf.and.bl_out(ib)%nrb.lt.nrafmax) THEN
          im(ib)=1
          bl_out(ib)%nrb=bl_out(ib)%nrb+1
          DO i=1,6
              IF (bl_out(ib)%cl(i).gt.0)im(bl_out(ib)%cl(i))=1
          END DO
      END IF
      IF (dabs(critere(ib)).lt.crit_deraf.and.bl_out(ib)%nrb.gt.0) THEN
          im(ib)=-1
          bl_out(ib)%nrb=bl_out(ib)%nrb-1
          DO i=1,6
              IF (bl_out(ib)%cl(i).gt.0)im(bl_out(ib)%cl(i))=1
          END DO
      END IF
  END DO
! on verifie que 2 blocs voisins n aient pas une difference d indice de raffinement >2      
 10  itest=0
  DO ib=1,nb_bloc
      DO i=1,6
          bcl=bl_out(ib)%cl(i)
          if(bcl.lt.-100000)bcl=-bcl-100000
          IF (bcl.gt.0) THEN
              IF ((bl_out(ib)%nrb-bl_out(bl_out(ib)%cl(i))%nrb).ge.2) THEN
                  itest=1
                  bl_out(bcl)%nrb=bl_out(bcl)%nrb+1
                  im(bcl)=1
                  idv=bcl
                  DO j=1,6
                      bclj=bl_out(idv)%cl(j)
                      if(bclj.lt.-100000)bclj=-bclj-100000
                      IF (bclj.gt.0) im(bclj)=1
                  END DO
              END IF
          END IF
      END DO
  END DO        
  IF(itest.ne.0) goto 10
        
! liste des voisins: -1/ le voisin est + grossier
!                     0/ le voisin est identique
!                     1/ le voisin est + fin
  DO ib=1,nb_bloc
      bl_out(ib)%nrf=0
      DO i=1,6
          bcl=bl_out(ib)%cl(i)
          if(bcl.lt.-100000)bcl=-bcl-100000
          IF (bcl.gt.0) THEN
              IF (bl_out(bcl)%nrb.gt.bl_out(ib)%nrb)bl_out(ib)%nrf(i)=1
              IF (bl_out(bcl)%nrb.lt.bl_out(ib)%nrb)bl_out(ib)%nrf(i)=-1
          END IF
      END DO
  END DO
         
! On refait la repartition des blocs dans chaque domaine
  DO ib=1,nb_bloc
       bl_out(ib)%idom=0
  END DO
  write( *,1011)nb_bloc,nbdom
 1011 format('Aggregation of ',i6,' blocks into',i3,' sub-domains')
  CALL splitdom_Z(bl_out,nb_bloc,nbdom, cb%msh%ndim) !!!!!!! DECOUPAGE PAR Z-ORDER

! Remaillage eventuel des blocs et projection des donnees
! Determination de la taille des tableaux du bloc
  DO ib=1,nb_bloc
      bl_out(ib)%dom=0
      DO i=1,6
          bcl=bl_out(ib)%cl(i)
          if(bcl.lt.-100000)bcl=-bcl-100000
          IF (bcl.gt.0) bl_out(ib)%dom(i)=bl_out(bcl)%idom
      END DO
      imac=1
      imav=1
      imaf=1
      DO i=1,6 ! ima: indice de maillage de l interface entre 2 blocs voisins
          bcl=bl_out(ib)%cl(i)
          if(bcl.lt.-100000)bcl=-bcl-100000
          IF (bl_out(ib)%idom.eq.bl_out(ib)%dom(i).and.ib.lt.bcl)imaf(i)=0
          IF (bl_out(ib)%idom.eq.bl_out(ib)%dom(i).and.ib.lt.bcl)imav(i)=0
          IF (bl_out(ib)%idom.eq.bl_out(ib)%dom(i).or.bcl.le.0)imac(i)=0
      END DO
      
      SELECT CASE (cb%msh%ndim)
        CASE (1) ! Cas 1D
          nx=bl_out(ib)%nx*2**(bl_out(ib)%nrb)
          ny=1
          nz=1
          bl_out(ib)%nb_cell   = nx+imac(1)+imac(2)
          bl_out(ib)%nb_vertex = (nx+1)*(ny+1)*2 
          bl_out(ib)%nb_face   = (nx-1)+nx*2*2 + imaf(1)+imaf(2) 

        CASE (2) ! Cas 2D
          nx=bl_out(ib)%nx*2**(bl_out(ib)%nrb)
          ny=bl_out(ib)%ny*2**(bl_out(ib)%nrb)
          nz=1
          bl_out(ib)%nb_cell = nx*ny + int(ny*2.**bl_out(ib)%nrf(1)*imac(1)) &
                                     + int(ny*2.**bl_out(ib)%nrf(2)*imac(2)) &
                                     + int(nx*2.**bl_out(ib)%nrf(3)*imac(3)) &
                                     + int(nx*2.**bl_out(ib)%nrf(4)*imac(4))
          bl_out(ib)%nb_vertex = (nx+1)*(ny+1)*2 &
                                     + bl_out(ib)%nrf(1)*(bl_out(ib)%nrf(1)+1)*ny*imav(1) &
                                     + bl_out(ib)%nrf(2)*(bl_out(ib)%nrf(2)+1)*ny*imav(2) &
                                     + bl_out(ib)%nrf(3)*(bl_out(ib)%nrf(3)+1)*nx*imav(3) &
                                     + bl_out(ib)%nrf(4)*(bl_out(ib)%nrf(4)+1)*nx*imav(4)
          bl_out(ib)%nb_face   = (nx-1)*ny+nx*(ny-1)+nx*ny*2 &
                                     + ny*imaf(1)+ny*imaf(2)+nx*imaf(3)+nx*imaf(4)  &
                                     + bl_out(ib)%nrf(1)*(bl_out(ib)%nrf(1)+1)*ny/2*imaf(1) &
                                     + bl_out(ib)%nrf(2)*(bl_out(ib)%nrf(2)+1)*ny/2*imaf(2) &
                                     + bl_out(ib)%nrf(3)*(bl_out(ib)%nrf(3)+1)*nx/2*imaf(3) &
                                     + bl_out(ib)%nrf(4)*(bl_out(ib)%nrf(4)+1)*nx/2*imaf(4)        
        CASE (3) ! Cas 3D
          nx=bl_out(ib)%nx*2**(bl_out(ib)%nrb)
          ny=bl_out(ib)%ny*2**(bl_out(ib)%nrb)
          nz=bl_out(ib)%nz*2**(bl_out(ib)%nrb)
          bl_out(ib)%nb_cell   = nx*ny*nz &
                                     + int(ny*nz*4.**bl_out(ib)%nrf(1)*imac(1)) &
                                     + int(ny*nz*4.**bl_out(ib)%nrf(2)*imac(2)) &
                                     + int(nx*nz*4.**bl_out(ib)%nrf(3)*imac(3)) &
                                     + int(nx*nz*4.**bl_out(ib)%nrf(4)*imac(4)) &
                                     + int(nx*ny*4.**bl_out(ib)%nrf(5)*imac(5)) &
                                     + int(nx*ny*4.**bl_out(ib)%nrf(6)*imac(6)) 
          bl_out(ib)%nb_vertex = (nx+1)*(ny+1)*(nz+1) &
                                     + bl_out(ib)%nrf(1)*(bl_out(ib)%nrf(1)+1)/2*(3*ny*nz+ny+nz)*imav(1) &
                                     + bl_out(ib)%nrf(2)*(bl_out(ib)%nrf(2)+1)/2*(3*ny*nz+ny+nz)*imav(2) &
                                     + bl_out(ib)%nrf(3)*(bl_out(ib)%nrf(3)+1)/2*(3*nx*nz+nx+nz)*imav(3) &
                                     + bl_out(ib)%nrf(4)*(bl_out(ib)%nrf(4)+1)/2*(3*nx*nz+nx+nz)*imav(4) &
                                     + bl_out(ib)%nrf(5)*(bl_out(ib)%nrf(5)+1)/2*(3*nx*ny+nx+ny)*imav(5) &
                                     + bl_out(ib)%nrf(6)*(bl_out(ib)%nrf(6)+1)/2*(3*nx*ny+nx+ny)*imav(6) 
          bl_out(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) &
                                     + bl_out(ib)%nrf(1)*(bl_out(ib)%nrf(1)+1)/2*3*ny*nz*imaf(1) &
                                     + bl_out(ib)%nrf(2)*(bl_out(ib)%nrf(2)+1)/2*3*ny*nz*imaf(2) &
                                     + bl_out(ib)%nrf(3)*(bl_out(ib)%nrf(3)+1)/2*3*nx*nz*imaf(3) &
                                     + bl_out(ib)%nrf(4)*(bl_out(ib)%nrf(4)+1)/2*3*nx*nz*imaf(4) &  
                                     + bl_out(ib)%nrf(5)*(bl_out(ib)%nrf(5)+1)/2*3*nx*ny*imaf(5) &  
                                     + bl_out(ib)%nrf(6)*(bl_out(ib)%nrf(6)+1)/2*3*nx*ny*imaf(6)   
                                                
      END SELECT
  END DO

! Creation de la liste des interfaces
  allocate(interdom(nbdom*(nbdom-1)/2))
  allocate(nbl(nbdom*(nbdom-1)/2))
  CALL inter(bl_out,nb_bloc,nbdom,interdom,nbl,cb%msh%ndim)

! Determination de petit
  petit=1.d-4
  DO ib=1,nb_bloc
      IF (hpetit/(bl_out(ib)%nx*2**bl_out(ib)%nrb).lt.petit) petit=hpetit/(bl_out(ib)%nx*2**bl_out(ib)%nrb)
      IF (hpetit/(bl_out(ib)%ny*2**bl_out(ib)%nrb).lt.petit.and.cb%msh%ndim.ne.1)&
      petit=hpetit/(bl_out(ib)%ny*2**bl_out(ib)%nrb)
      IF (hpetit/(bl_out(ib)%nz*2**bl_out(ib)%nrb).lt.petit.and.cb%msh%ndim.eq.3)&
      petit=hpetit/(bl_out(ib)%nz*2**bl_out(ib)%nrb)
  END DO
  petit=petit*0.1d0
  write( *,1818)petit,hpetit
  1818     format('CERF_AMR-INFO Small dimensions :',2e13.5)

! Boucle sur les domaines
  DO idom=1,nbdom
      ! Initialisations
      calc_out%msh%ndim  =cb%msh%ndim
      calc_out%msh%nb_dom=nbdom
      nullify(calc_out%msh%list_vertex)
      nullify(calc_out%msh%list_face)
      nullify(calc_out%msh%list_cell)
      nullify(calc_out%msh%list_send)
      nullify(calc_out%msh%list_recv)
      CALL empty_mesh(calc_out%msh,calc_out%msh%nb_dom,idom)
      CALL copy_sol(cb,calc_out)
      IF (associated(calc_out%list_bloc)) THEN
          deallocate(calc_out%list_bloc)
          nullify(calc_out%list_bloc)
      END IF
      nullify(calc_out%list_bloc)

      ! Determination de la taille des tableaux du domaine
      calc_out%msh%nb_cell  = 0
      calc_out%msh%nb_vertex= 0
      calc_out%msh%nb_face  = 0
      calc_out%nb_bloc=0
         
      DO ib=1,nb_bloc
          IF (bl_out(ib)%idom.eq.idom) THEN
              calc_out%nb_bloc=calc_out%nb_bloc+1
              bl_out(ib)%first_cell=calc_out%msh%nb_cell
              calc_out%msh%nb_cell  = calc_out%msh%nb_cell  +bl_out(ib)%nb_cell
              bl_out(ib)%first_face=calc_out%msh%nb_face
              calc_out%msh%nb_face  = calc_out%msh%nb_face  +bl_out(ib)%nb_face
              bl_out(ib)%first_vertex=calc_out%msh%nb_vertex
              calc_out%msh%nb_vertex= calc_out%msh%nb_vertex+bl_out(ib)%nb_vertex
          END IF
      END DO
         
      allocate(calc_out%list_bloc(calc_out%nb_bloc))
      allocate(calc_out%msh%list_cell(calc_out%msh%nb_cell))
      allocate(calc_out%msh%list_vertex(calc_out%msh%nb_vertex))
      allocate(calc_out%msh%list_face(calc_out%msh%nb_face))
      allocate(calc_out%msh%list_send(calc_out%msh%nb_dom))
      allocate(calc_out%msh%list_recv(calc_out%msh%nb_dom))
      DO i=1,idom-1
          k=idom-i+(i-1)*(2*nbdom-i)/2 
          IF (nbl(k).ne.0) allocate(calc_out%msh%list_recv(i)%L(interdom(k)%L(4*nbl(k)-1)))
          IF (nbl(k).ne.0) allocate(calc_out%msh%list_send(i)%L(interdom(k)%L(4*nbl(k)-0)))
      END DO
      DO i=idom+1,nbdom
          k=i-idom+(idom-1)*(2*nbdom-idom)/2
          IF (nbl(k).ne.0) allocate(calc_out%msh%list_recv(i)%L(interdom(k)%L(4*nbl(k)-0)))
          IF (nbl(k).ne.0) allocate(calc_out%msh%list_send(i)%L(interdom(k)%L(4*nbl(k)-1)))
      END DO

! Maillage des blocs du domaine
      icb=0
      DO ib=1,nb_bloc
          IF (bl_out(ib)%idom.eq.idom) THEN
              icb=icb+1
              calc_out%list_bloc(icb)=bl_out(ib)
              DO i=1,6
                  blv(i)=bl_null
                  bcl=bl_out(ib)%cl(i)
                  if(bcl.lt.-100000)bcl=-bcl-100000
                  IF (bcl.gt.0)blv(i)=bl_out(bcl)
              END DO
              SELECT CASE (cb%msh%ndim)
                 CASE (1)
                   CALL mesh1d(bl_out(ib),calc_out%msh,ib,nbdom,blv,interdom)
                 CASE (2)
                  CALL mesh2d(bl_out(ib),calc_out%msh,ib,nbdom,blv,interdom)
                 CASE (3)
                  CALL mesh3d(bl_out(ib),calc_out%msh,ib,nbdom,blv,interdom)
              END SELECT
          END IF
      END DO ! Fin boucle sur les blocs
      CALL init_mesh(calc_out%msh) 
      
! On verifie les orientations de face
      iv=0
      DO i=1,calc_out%msh%nb_face
          ic1=calc_out%msh%list_face(i)%ic1
          d=(calc_out%msh%list_face(i)%center%x-calc_out%msh%list_cell(ic1)%center%x)*&
             calc_out%msh%list_face(i)%vnorm%v(1)+&
            (calc_out%msh%list_face(i)%center%y-calc_out%msh%list_cell(ic1)%center%y)*&
             calc_out%msh%list_face(i)%vnorm%v(2)+&
            (calc_out%msh%list_face(i)%center%z-calc_out%msh%list_cell(ic1)%center%z)*&
            calc_out%msh%list_face(i)%vnorm%v(3)         
          IF (d.le.0.) THEN
              calc_out%msh%list_face(i)%vnorm%v=-calc_out%msh%list_face(i)%vnorm%v
              itemp=calc_out%msh%list_face(i)%vertex(2)
              calc_out%msh%list_face(i)%vertex(2)=calc_out%msh%list_face(i)%vertex(calc_out%msh%list_face(i)%nbvertex)
              calc_out%msh%list_face(i)%vertex(calc_out%msh%list_face(i)%nbvertex)=itemp
              iv=iv+1
          END IF
      END DO
      IF (iv.ne.0) CALL init_mesh(calc_out%msh)

! On elimine les sommets surabondants
      nvi=calc_out%msh%nb_vertex
      CALL trivertex(calc_out%msh)
      write( *,1111)idom,nvi-calc_out%msh%nb_vertex,nvi,calc_out%msh%nb_vertex
 1111 format('CERF_AMR-INFO Domain',i3,' elimination of ',i7,' vertexes out of',i7,', i.e.',i7,' vertexes.')
        
! Initialisation sur tout le domaine
      nbtotcello=0
      nbtotcinto=0
      CALL empty_mesh(calc_ini%msh,nbdom,idom)
      CALL empty_sol(calc_ini%list_sol)
      CALL get_calc(calc_ini)
     
! Copie de l objet calcul              
      calc_out%tmin           = calc_ini%tmin
      calc_out%tmax           = calc_ini%tmax
      calc_out%dt             = calc_ini%dt
      calc_out%pdtp           = calc_ini%pdtp
      calc_out%cfl            = calc_ini%cfl
      calc_out%vmax           = calc_ini%vmax
      calc_out%isauv          = calc_ini%isauv
      calc_out%istart         = calc_ini%istart
      calc_out%ipas           = calc_ini%ipas
      calc_out%iordre_e       = calc_ini%iordre_e
      calc_out%iordre_t       = calc_ini%iordre_t
      calc_out%imeth_int_temps= calc_ini%imeth_int_temps
      calc_out%point_sondes   = calc_ini%point_sondes
      calc_out%nb_sondes      = calc_ini%nb_sondes
      calc_out%num_point_min  = calc_ini%num_point_min
      calc_out%ityp_flux      = calc_ini%ityp_flux
      calc_out%ilimiteur      = calc_ini%ilimiteur
      calc_out%nbds           = calc_ini%nbds
      calc_out%nrma           = calc_ini%nrma
      calc_out%vcde           = calc_ini%vcde
      calc_out%vcra           = calc_ini%vcra
      calc_out%nb_sol         = calc_ini%nb_sol
      calc_out%list_sol       = calc_ini%list_sol
            
           
      DO ib=1,calc_out%nb_bloc
          ibloc=calc_out%list_bloc(ib)%numbl
          CALL proddl(list_w_bloc(ibloc),bl_ini(ibloc),bl_out(ibloc),calc_out,cb%msh%ndim)
      END DO
!---------------------------
! Si presence de solides
      IF (calc_out%nb_sol .NE. 0) THEN
          DO i = 1,calc_out%msh%nb_cell
            calc_out%msh%list_cell(i)%w%chi=0.d0
          enddo
          DO j = 1,calc_out%nb_sol
               k=0
               DO i = 1,calc_out%msh%nb_cell
                    IF (          &
                    calc_out%msh%list_cell(i)%center%x<calc_out%list_sol(j)%xmax_s  &
                    .AND. calc_out%msh%list_cell(i)%center%x>calc_out%list_sol(j)%xmin_s  &
                    .AND. calc_out%msh%list_cell(i)%center%y<calc_out%list_sol(j)%ymax_s  &
                    .AND. calc_out%msh%list_cell(i)%center%y>calc_out%list_sol(j)%ymin_s  &
                    .AND. calc_out%msh%list_cell(i)%center%z<calc_out%list_sol(j)%zmax_s  &
                    .AND. calc_out%msh%list_cell(i)%center%z>calc_out%list_sol(j)%zmin_s) &
                    CALL raycasting(calc_out%list_sol(j),calc_out%msh%list_cell(i)%center%x,&
                                                     calc_out%msh%list_cell(i)%center%y,&
                                                     calc_out%msh%list_cell(i)%center%z,&
                                                     calc_out%msh%list_cell(i)%w%chi)
                    IF (calc_out%msh%list_cell(i)%w%chi .EQ. calc_out%list_sol(j)%id) k=k+1
               END DO
               calc_out%list_sol(j)%nb_cell = k
          END DO

          DO j = 1,calc_out%nb_sol
               k=0
               DO i = 1,calc_out%msh%nb_cell
                    IF (   calc_out%msh%list_cell(i)%idom .EQ. calc_out%msh%numdom         &
                    .AND.  calc_out%msh%list_cell(i)%center%x<calc_out%list_sol(j)%xmax_s  &
                    .AND.  calc_out%msh%list_cell(i)%center%x>calc_out%list_sol(j)%xmin_s  &
                    .AND.  calc_out%msh%list_cell(i)%center%y<calc_out%list_sol(j)%ymax_s  &
                    .AND.  calc_out%msh%list_cell(i)%center%y>calc_out%list_sol(j)%ymin_s  &
                    .AND.  calc_out%msh%list_cell(i)%center%z<calc_out%list_sol(j)%zmax_s  &
                    .AND.  calc_out%msh%list_cell(i)%center%z>calc_out%list_sol(j)%zmin_s  &
                    .AND. (calc_out%msh%list_cell(i)%w%chi .EQ. calc_out%list_sol(j)%id)) THEN
                         k=k+1
                         !calc_out%list_sol(j)%list_cells(k) = i
                         calc_out%msh%list_cell(i)%w%chi = calc_out%list_sol(j)%id
                         calc_out%msh%list_cell(i)%w%vprim%v(1)=calc_out%list_sol(j)%rho
                         CALL prim2bal(calc_out%msh%list_cell(i)%w)
                    END IF
               END DO
          END DO
          if(idom.eq.1)call save_obj(calc_out)
      END IF     
      

! Sauvegarde
      calc_out%msh%numdom=idom
      DO i=1,calc_out%msh%nb_cell
          calc_out%msh%list_cell(i)%w%vbal%v=0.d0
          calc_out%msh%list_cell(i)%balance%v=0.d0
          calc_out%msh%list_cell(i)%fp%v=0.d0
          calc_out%msh%list_cell(i)%w%gprim(1)%v=0.d0
          calc_out%msh%list_cell(i)%w%gprim(2)%v=0.d0
          calc_out%msh%list_cell(i)%w%gprim(3)%v=0.d0
          calc_out%msh%list_cell(i)%w%enti=0.d0
          calc_out%msh%list_cell(i)%w%entf=0.d0
          calc_out%msh%list_cell(i)%w%pe  =0.d0
          CALL prim2bal(calc_out%msh%list_cell(i)%w)
      END DO

      delta  =calc_out%tmax-calc_out%tmin
      calc_out%tmin  =calc_out%tmin-delta
      calc_out%tmax  =calc_out%tmax-delta     
 
      CALL save_calc(calc_out)
      CALL empty_mesh(calc_out%msh,0,0)  ! liberation de la mémoire
      CALL empty_mesh(calc_ini%msh,0,0)  ! liberation de la mémoire
      CALL empty_sol(calc_ini%list_sol)
      CALL empty_sol(calc_out%list_sol)
    
  END DO
  
  
  nbtotcinti=0
  nbtotcinto=0
  nbtotcelli=0   
  nbtotcello=0   
  DO ib=1,nb_bloc
      SELECT CASE (cb%msh%ndim)
        CASE (1)
          nb_cell_bloc=bl_ini(ib)%nx*2**bl_ini(ib)%nrb
        CASE (2)
          nb_cell_bloc=bl_ini(ib)%nx*2**bl_ini(ib)%nrb*bl_ini(ib)%ny*2**bl_ini(ib)%nrb
        CASE (3)
          nb_cell_bloc=bl_ini(ib)%nx*2**bl_ini(ib)%nrb* &
                       bl_ini(ib)%ny*2**bl_ini(ib)%nrb* &
                       bl_ini(ib)%nz*2**bl_ini(ib)%nrb
      END SELECT
      nbtotcinti=nbtotcinti+nb_cell_bloc
      nbtotcelli=nbtotcinti+bl_ini(ib)%nb_cell
      SELECT CASE (cb%msh%ndim)
        CASE (1)
          nb_cell_bloc=bl_out(ib)%nx*2**bl_out(ib)%nrb
        CASE (2)
          nb_cell_bloc=bl_out(ib)%nx*2**bl_out(ib)%nrb* &
                       bl_out(ib)%ny*2**bl_out(ib)%nrb
        CASE (3)
         nb_cell_bloc=bl_out(ib)%nx*2**bl_out(ib)%nrb* &
                      bl_out(ib)%ny*2**bl_out(ib)%nrb* &
                      bl_out(ib)%nz*2**bl_out(ib)%nrb
      END SELECT
      nbtotcinto=nbtotcinto+nb_cell_bloc
      nbtotcello=nbtotcinto+bl_ini(ib)%nb_cell
  END DO
  write(*,*)'Total number of cells,    before', nbtotcelli, '  after', nbtotcello
  write(*,*)'Number of interior cells, before', nbtotcinti, '  after', nbtotcinto
!===========================   FIN DE LA ROUTINE    ====================

  END PROGRAM cerf_amr
!DOCFIN
