TABLE OF CONTENTS


UTI/cerf_amr [ Modules ]

[ Top ] [ Modules ]

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