cerf_amr.f90

### TABLE OF CONTENTS

- 1. UTI/cerf_amr

## 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