mesh_interpolation.f90

Go to the documentation of this file.

MODULE mesh_interpolation

  PUBLIC :: mesh_interpol
  PRIVATE
CONTAINS
  SUBROUTINE mesh_interpol

    USE chaine_caractere
    USE def_type_mesh
    USE my_util
    USE create_comm
    USE interpolation_tools
    USE metis_sfemans
    USE prep_maill
    USE sub_plot
    USE restart

    IMPLICIT NONE

    TYPE(mesh_type)                       :: p1_mesh_glob, p2_mesh_glob, p2_c0_mesh_glob_temp

    TYPE(mesh_type), TARGET               :: h_mesh, phi_mesh
    TYPE(mesh_type), TARGET               :: h_mesh_glob, phi_mesh_glob
    TYPE(mesh_type), POINTER              :: h_mesh_in, h_mesh_out, phi_mesh_in, phi_mesh_out

    TYPE(mesh_type), TARGET               :: vv_mesh, pp_mesh
    TYPE(mesh_type), TARGET               :: vv_mesh_glob, pp_mesh_glob
    TYPE(mesh_type), POINTER              :: vv_mesh_in, vv_mesh_out, pp_mesh_in, pp_mesh_out

    TYPE(mesh_type), TARGET               :: temp_mesh
    TYPE(mesh_type), TARGET               :: temp_mesh_glob
    TYPE(mesh_type), POINTER              :: temp_mesh_in, temp_mesh_out

    CHARACTER(len=200)                    :: old_filename, old_directory
    CHARACTER(len=200)                    :: new_filename, new_directory
    CHARACTER(len=200)                    :: file_name, directory
    CHARACTER(len=200)                    :: file_name_m, directory_m
    LOGICAL                               :: if_momentum, if_mass, if_induction, if_energy, &
         old_is_form, new_is_form, iformatted, is_form_m, mono_in, mono_out, is_in, &
         dom_h_larger_dom_ns, dom_temp_larger_dom_ns, &
         if_temperature_in, if_temperature_out, if_temperature, if_level_set, inter_mesh, &
         rw_ns, rw_temp, rw_mxw, check_plt, test, if_select

    INTEGER                               :: nb_s, nb_f, nb_procs, rank, type_fe_h, type_fe_phi, &
         nb_dom_ns, nb_dom_temp, nb_dom_h, nb_dom_phi, nsize, code, m, &
         m_max_c, rank_s, nb_inter, nb_inter_mu, nb_inter_v_t, &
         k, kp, i, nb_mode, rang_ns_s, rang_temp_s, rang_s, l, lblank, nb_fluid
    INTEGER                               :: nb_fic, index_start
    REAL(KIND=8)                          :: time_h, time_u, time_t, max_vel
    TYPE(periodic_data)                   :: my_periodic

    INTEGER, DIMENSION(:), ALLOCATABLE    :: list_inter, list_inter_temp, list_inter_v_t, &
         list_inter_mu, list_inter_h_phi, list_dom_h, & 
         list_dom_ns, list_dom_temp, list_dom_temp_in, list_dom_phi, list_dom_h_in, list_dom, part, &
         list_mode, controle_h, controle_phi, &
         controle_vv, controle_pp, controle_temp, temp_in_to_new, h_in_to_new, &
         l_t_g_vv, l_t_g_pp, l_t_g_h, l_t_g_phi, l_t_g_temp

    REAL(KIND=8), DIMENSION(:,:,:),   POINTER     :: hn, hn1, bn, bn1, phin, phin1
    REAL(KIND=8), DIMENSION(:,:,:),   POINTER     :: hn_glob, hn1_glob, bn_glob, bn1_glob, phin_glob, phin1_glob
    REAL(KIND=8), DIMENSION(:,:,:),   POINTER     :: un, un_m1, pn, pn_m1
    REAL(KIND=8), DIMENSION(:,:,:),   POINTER     :: incpn, incpn_m1, tempn, tempn_m1
    REAL(KIND=8), DIMENSION(:,:,:,:), POINTER     :: level_setn, level_setn_m1
    REAL(KIND=8), DIMENSION(:,:,:),   POINTER     :: un_glob, un_m1_glob, pn_glob, pn_m1_glob
    REAL(KIND=8), DIMENSION(:,:,:),   POINTER     :: incpn_glob, incpn_m1_glob, tempn_glob, tempn_m1_glob
    REAL(KIND=8), DIMENSION(:,:,:,:), POINTER     :: level_setn_glob, level_setn_m1_glob

    REAL(KIND=8), DIMENSION(:,:,:),   POINTER     :: hn_in, hn1_in, bn_in, bn1_in, phin_in, phin1_in
    REAL(KIND=8), DIMENSION(:,:,:),   POINTER     :: hn_out, hn1_out, bn_out, bn1_out, phin_out, phin1_out
    REAL(KIND=8), DIMENSION(:,:,:),   POINTER     :: un_in, un_m1_in, pn_in, pn_m1_in
    REAL(KIND=8), DIMENSION(:,:,:),   POINTER     :: incpn_in, incpn_m1_in, tempn_in, tempn_m1_in
    REAL(KIND=8), DIMENSION(:,:,:,:), POINTER     :: level_setn_in, level_setn_m1_in
    REAL(KIND=8), DIMENSION(:,:,:),   POINTER     :: un_out, un_m1_out, pn_out, pn_m1_out
    REAL(KIND=8), DIMENSION(:,:,:),   POINTER     :: incpn_out, incpn_m1_out, tempn_out, tempn_m1_out
    REAL(KIND=8), DIMENSION(:,:,:,:), POINTER     :: level_setn_out, level_setn_m1_out

    CHARACTER(len=3)             :: type_pb, tit_m, tit_s, tit
    CHARACTER(len=4)             :: tit_part
    CHARACTER(len=200)           :: mesh_part_name
    INTEGER                      :: petsc_rank
    LOGICAL                      :: if_read_partition

#include "petsc/finclude/petsc.h"
    mpi_comm                                        :: comm_cart
    mpi_comm, DIMENSION(:), POINTER                 :: comm_one_d, comm_one_d_ns, comm_one_d_temp, coord_cart
    petscerrorcode                                  :: ierr

    CALL petscinitialize(petsc_null_character,ierr)
    CALL mpi_comm_rank(petsc_comm_world, petsc_rank, code)

    OPEN(unit=22, file='data_interpol', status='unknown', access='sequential')

    !===Compute is_in (is_in=.TRUE. if Petsc -> NonPetsc; is_in=.FALSE. if NonPetsc -> Petsc)
    CALL which_pb(is_in)

    IF (is_in) THEN
       CALL read_until(22, '===Number of processors in meridian section (Input)')
       READ(22,*) nb_s
       CALL read_until(22, '===Problem type (Input): (nst, mxw, mhd, fhd)')
       READ(22,*) type_pb
       CALL find_string(22, '===Is there an input temperature field?', test)
       IF (test) THEN
          READ (22, *) if_temperature_in
       ELSE
          if_temperature_in = .false.
       END IF
       if_temperature = if_temperature_in
       inter_mesh = .false.
       if_read_partition = .true.
    ELSE
       CALL read_until(22, '===Number of processors in meridian section (Output)')
       READ(22,*) nb_s
       CALL read_until(22, '===Problem type (Output): (nst, mxw, mhd, fhd)')
       READ(22,*) type_pb
       CALL find_string(22, '===Is there an output temperature field?', test)
       IF (test) THEN
          READ (22, *) if_temperature_out
       ELSE
          if_temperature_out = .false.
       END IF
       if_temperature = if_temperature_out
       CALL read_until(22, '===Should data be interpolated on new mesh? (True/False)')
       READ(22,*) inter_mesh
       if_read_partition = .false.
    END IF

    !===Set rw_ns, rw_mxw and rw_temp (function of Problem type (Input))
    CALL read_until(22, '===Problem type (Input): (nst, mxw, mhd, fhd)')
    READ(22,*) tit
    IF (tit=='nst') THEN
       rw_ns = .true.
       rw_mxw = .false.
    ELSE IF ((tit=='mhd').OR.(tit=='fhd')) THEN
       rw_ns = .true.
       rw_mxw = .true.
    ELSE IF (tit=='mxw') THEN
       rw_ns = .false.
       rw_mxw = .true.
    ELSE
       CALL error_petsc('BUG in interpol, type_pb (Input) not correct')
    END IF
    CALL find_string(22, '===Is there an input temperature field?', test)
    IF (test) THEN
       READ (22, *) rw_temp
    ELSE
       rw_temp = .false.
    END IF

    !===Default parameters
    nb_f = 1

    !===Check construction with plotmtv
    CALL find_string(22, '===Check construction with plotmtv? (True/False)', test)
    IF (test) THEN
       READ (22, *) check_plt
    ELSE
       check_plt  = .false.
    END IF

    !===Number of time steps to process
    CALL read_until(22, '===How many files (times steps) should be converted?')
    READ(22, *) nb_fic

    !===Starting index
    CALL find_string(22, '===What is the starting index? (suite*_I*index.mesh*), index is an integer in [1,999]', test)
    IF (test) THEN
       READ (22, *) index_start
    ELSE
       index_start  = 1
    END IF

    !===Number of Fourier modes
    CALL read_until(22, '===Number of Fourier modes')
    READ(22, *) nb_mode
    ALLOCATE(list_mode(nb_mode))
    CALL find_string(22, '===Select Fourier modes? (true/false)',test)
    IF (test) THEN
       READ(22,*) if_select
       IF (if_select) THEN
          CALL read_until(22, '===List of Fourier modes (if select_mode=.TRUE.)')
          READ(22,*) list_mode
       ELSE
          list_mode = (/ (m, m=0, nb_mode-1) /)
       END IF
    ELSE
       list_mode = (/ (m, m=0, nb_mode-1) /)
    END IF

    !===Input and Output mesh files
    CALL read_until(22, '===Directory and name of input mesh file') 
    READ(22,*) old_directory, old_filename
    CALL read_until(22, '===Is input mesh file formatted (true/false)?')
    READ(22,*) old_is_form
    IF (is_in) THEN !Interpolation is done at second step only
       new_directory = old_directory
       new_filename = old_filename
       new_is_form = old_is_form
    ELSE
       CALL read_until(22, '===Directory and name of output mesh file') 
       READ(22,*) new_directory, new_filename
       CALL read_until(22, '===Is output mesh file formatted (true/false)?')
       READ(22,*) new_is_form
    END IF

    !===Data for NS
    IF ( (type_pb == 'nst') .OR. (type_pb == 'mhd') .OR. (type_pb == 'fhd')) THEN
       CALL read_until(22, '===Number of subdomains in Navier-Stokes mesh')
       READ(22,*) nb_dom_ns 
       ALLOCATE(list_dom_ns(nb_dom_ns))
       CALL read_until(22, '===List of subdomains for Navier-Stokes mesh')
       READ (22, *)  list_dom_ns
       CALL find_string(22, '===Is there a level set?', test)
       IF (test) THEN
          READ (22, *) if_level_set
          IF (if_level_set) THEN
             CALL read_until(22, '===How many fluids?')
             READ(22,*) nb_fluid
          END IF
       ELSE
          if_level_set = .false.
       END IF
       IF (type_pb == 'nst') ALLOCATE(list_dom_h(0), list_dom_phi(0))
    END IF

    !===Data for temperature
    IF (if_temperature) THEN
       CALL read_until(22, '===Number of subdomains in temperature mesh')
       READ(22,*) nb_dom_temp 
       ALLOCATE(list_dom_temp_in(nb_dom_temp), list_dom_temp(nb_dom_temp), temp_in_to_new(nb_dom_temp))
       CALL read_until(22, '===List of subdomains for temperature mesh')
       READ (22, *)  list_dom_temp_in
       CALL find_string(21, '===Number of interfaces between velocity and temperature only domains (for nst applications)', test)
       IF (test) THEN
          READ(21,*) nb_inter_v_t
          ALLOCATE(list_inter_v_t(nb_inter_v_t))
          CALL read_until(21, '===List of interfaces between velocity and temperature only domains (for nst applications)')
          READ(21,*) list_inter_v_t
       ELSE
          ALLOCATE(list_inter_v_t(0))
       END IF
    ELSE
       ALLOCATE(list_dom_temp(0))
    END IF

    !===Data for Maxwell
    IF ( (type_pb == 'mhd') .OR. (type_pb == 'mxw') .OR. (type_pb == 'fhd')) THEN
       CALL read_until(22, '===Type of finite element for magnetic field')
       READ (22, *) type_fe_h
       CALL read_until(22, '===Number of subdomains in magnetic field (H) mesh')
       READ(22,*) nb_dom_h  ! number of sub_domains for H
       ALLOCATE(list_dom_h_in(nb_dom_h), list_dom_h(nb_dom_h), h_in_to_new(nb_dom_h)) ! JLG/AR Nov 17 2008
       CALL read_until(22, '===List of subdomains for magnetic field (H) mesh')
       READ(22,*) list_dom_h_in

       !==========Interfaces H/H==========================!
       CALL read_until(22, '===Number of interfaces in H mesh')
       READ(22,*) nb_inter_mu
       ALLOCATE(list_inter_mu(nb_inter_mu))
       IF (nb_inter_mu>0) THEN
          CALL read_until(22, '===List of interfaces in H mesh')
          READ(22,*) list_inter_mu
       END IF

       !==========Subdomains for phi======================!
       CALL find_string(22, '===Number of subdomains in magnetic potential (phi) mesh', test)
       IF (test) THEN
          READ (22, *) nb_dom_phi
       ELSE
          nb_dom_phi = 0
       END IF
       ALLOCATE(list_dom_phi(nb_dom_phi))
       IF (nb_dom_phi>0) THEN
          !==========List of subdomains for phi======================!
          CALL read_until(22, '===List of subdomains for magnetic potential (phi) mesh')
          READ(22,*) list_dom_phi

          !==========H/phi interface=========================!
          CALL read_until(22, '===Number of interfaces between H and phi')
          READ(22,*) nb_inter  ! number of interfaces between H and phi
          ALLOCATE(list_inter_h_phi(nb_inter))
          IF (nb_inter>0) THEN
             CALL read_until(22, '===List of interfaces between H and phi')
             READ(22,*) list_inter_h_phi
          END IF

          !==========Finite element type=====================!
          CALL read_until(22, '===Type of finite element for scalar potential')
          READ(22,*) type_fe_phi
       ELSE
          ALLOCATE(list_inter_h_phi(0))
          type_fe_phi = -1
       END IF
       IF (type_pb == 'mxw') THEN
          ALLOCATE(list_dom_ns(0))
       END IF
    END IF

    CALL find_string(22, '===How many pieces of periodic boundary?', test)
    IF (test) THEN
       READ (22, *) my_periodic%nb_periodic_pairs
    ELSE
       my_periodic%nb_periodic_pairs = 0
    END IF

    IF (my_periodic%nb_periodic_pairs.GE.1) THEN
       ALLOCATE(my_periodic%list_periodic(2,my_periodic%nb_periodic_pairs))
       ALLOCATE(my_periodic%vect_e(2,my_periodic%nb_periodic_pairs))
       CALL read_until(22, '===Indices of periodic boundaries and corresponding vectors')
       DO k = 1, my_periodic%nb_periodic_pairs
          READ(22,*) my_periodic%list_periodic(:,k), my_periodic%vect_e(:,k)
       END DO
    END IF
    CLOSE(22)

    !===Creation of logicals for equations
    if_mass = if_level_set
    if_momentum = type_pb=='nst' .OR. type_pb=='mhd' .OR. type_pb=='fhd'
    if_induction = type_pb=='mxw' .OR. type_pb=='mhd' .OR. type_pb=='fhd' ! No mxx: for mxx, put mhd
    if_energy = if_temperature

    !===Check that vv_mesh is a subset of H_mesh
    IF (if_induction) THEN
       IF (if_momentum) THEN
          IF (SIZE(list_dom_h) < SIZE(list_dom_ns)) THEN
             WRITE(*,*) ' BUG: NS must be a subset of Maxwell ' 
             stop
          END IF
          DO k = 1, nb_dom_ns
             ! JLG/AR Nov 17 2008 
             IF (minval(abs(list_dom_h_in - list_dom_ns(k))) /= 0) THEN
                WRITE(*,*) ' BUG : NS must be a subset of Maxwell '
                stop
             END IF
             DO kp = 1, nb_dom_h
                IF (list_dom_h_in(kp) == list_dom_ns(k)) EXIT  
             END DO
             h_in_to_new(k) = kp
             ! JLG/AR Nov 17 2008 
             list_dom_h(k) = list_dom_ns(k)
          END DO
          m = nb_dom_ns 
          DO k = 1, nb_dom_h
             IF (minval(abs(list_dom_h_in(k) - list_dom_ns)) == 0) cycle 
             m = m + 1
             ! JLG/AR Nov 17 2008 
             h_in_to_new(m) = k
             ! JLG/AR Nov 17 2008 
             list_dom_h(m) = list_dom_h_in(k)
          END DO
          IF (m/=nb_dom_h) THEN
             WRITE(*,*) ' BUG : m/=nb_dom_H ' 
             stop
          END IF
          IF (nb_dom_h > nb_dom_ns) THEN
             dom_h_larger_dom_ns = .true.
          ELSE
             dom_h_larger_dom_ns = .false.
          END IF
       ELSE
          ! JLG/AR Nov 17 2008 
          DO k = 1, nb_dom_h
             h_in_to_new(k) = k
          END DO
          ! JLG/AR Nov 17 2008 
          list_dom_h = list_dom_h_in
       END IF
    END IF

    !===Check that vv_mesh is a subset of temp_mesh
    IF (if_energy) THEN
       IF (SIZE(list_dom_temp) < SIZE(list_dom_ns)) THEN
          WRITE(*,*) ' BUG: NS must be a subset of temp ' 
          stop
       END IF
       DO k = 1, nb_dom_ns
          IF (minval(abs(list_dom_temp_in - list_dom_ns(k))) /= 0) THEN
             WRITE(*,*) ' BUG : NS must be a subset of temp '
             stop
          END IF
          DO kp = 1, nb_dom_temp
             IF (list_dom_temp_in(kp) == list_dom_ns(k)) EXIT  
          END DO
          temp_in_to_new(k) = kp
          list_dom_temp(k) = list_dom_ns(k)
       END DO
       m = nb_dom_ns 
       DO k = 1, nb_dom_temp
          IF (minval(abs(list_dom_temp_in(k) - list_dom_ns)) == 0) cycle 
          m = m + 1
          temp_in_to_new(m) = k
          list_dom_temp(m) = list_dom_temp_in(k)
       END DO
       IF (m/=nb_dom_temp) THEN
          WRITE(*,*) ' BUG : m/=nb_dom_temp ' 
          stop
       END IF
       IF (nb_dom_temp > nb_dom_ns) THEN
          dom_temp_larger_dom_ns = .true.
       ELSE
          dom_temp_larger_dom_ns = .false.
       END IF
    END IF

    !===Check that temp_mesh is a subset of H_mesh
    IF (if_induction .AND. if_energy) THEN
       IF (SIZE(list_dom_h) < SIZE(list_dom_temp)) THEN
          WRITE(*,*) ' BUG: temp must be a subset of Maxwell ' 
          stop
       END IF
       DO k = 1, nb_dom_temp
          IF (minval(abs(list_dom_h - list_dom_temp(k))) /= 0) THEN
             WRITE(*,*) ' BUG: temp must be a subset of Maxwell ' 
             stop
          END IF
       END DO
    END IF

    !===Create interfaces in meshes
    IF (if_momentum .AND. (.NOT. if_induction)) THEN
       IF (if_energy) THEN
          nsize = SIZE(list_dom_temp)
          ALLOCATE(list_dom(nsize))
          list_dom = list_dom_temp
          ALLOCATE(list_inter(SIZE(list_inter_v_t)))
          list_inter = list_inter_v_t
       ELSE
          nsize = SIZE(list_dom_ns)
          ALLOCATE(list_dom(nsize))
          list_dom = list_dom_ns
          ALLOCATE(list_inter(0))
       END IF
    ELSE
       nsize = SIZE(list_dom_h)+SIZE(list_dom_phi)
       ALLOCATE(list_dom(nsize))
       list_dom(1:SIZE(list_dom_h)) = list_dom_h
       IF (SIZE(list_dom_phi)>0) THEN
          list_dom(SIZE(list_dom_h)+1:) = list_dom_phi
       END IF
       nsize = SIZE(list_inter_mu)+SIZE(list_inter_h_phi)
       ALLOCATE(list_inter(nsize))
       IF (SIZE(list_inter_mu)>0) THEN
          list_inter(1:SIZE(list_inter_mu)) = list_inter_mu
       END IF
       IF (SIZE(list_inter_h_phi)>0) THEN 
          list_inter(SIZE(list_inter_mu)+1:) = list_inter_h_phi
       END IF
    END IF
    IF (if_energy) THEN
       ALLOCATE(list_inter_temp(0))
    END IF

    !===Directory, file name and format
    IF (is_in) THEN
       directory  = old_directory
       file_name  = old_filename
       iformatted = old_is_form
    ELSE
       directory  = new_directory
       file_name  = new_filename
       iformatted = new_is_form
    END IF

    IF (is_in) THEN
       directory_m = new_directory
       file_name_m = new_filename
       is_form_m   = new_is_form
    ELSE
       directory_m = old_directory
       file_name_m = old_filename
       is_form_m   = old_is_form
    END IF

    !===Check coherence
    CALL mpi_comm_size(petsc_comm_world, nb_procs, ierr)
    IF ( nb_s*nb_f /= nb_procs ) THEN
       CALL error_petsc('BUG: nb_S * nb_F /= nb_procs')
    END IF

    CALL create_cart_comm((/nb_s, nb_f/),comm_cart,comm_one_d,coord_cart)
    CALL mpi_comm_size(comm_one_d(2),nb_procs,code)
    CALL mpi_comm_rank(comm_one_d(2),rank,code)
    IF (nb_f/=nb_procs) THEN
       CALL error_petsc('BUG in INIT, nb_procs_F/=nb_procs')
    END IF

    IF ( (rw_ns) .AND. (type_pb=='mxw')) THEN
       CALL error_petsc('WARNING: type_pb and/or rw_ns might be wrong')
    END IF
    IF ( (rw_mxw) .AND. ((type_pb=='nst')) ) THEN
       CALL error_petsc('WARNING: type_pb and/or rw_mxw might be wrong')
    END IF

    !===Prepare meshes and pointers
    CALL load_dg_mesh_free_format(directory, file_name, list_dom, list_inter, 1, p1_mesh_glob, iformatted)
    CALL load_dg_mesh_free_format(directory, file_name, list_dom, list_inter, 2, p2_mesh_glob, iformatted)
    IF (if_energy) THEN
       CALL load_dg_mesh_free_format(directory, file_name, list_dom, &
            list_inter_temp, 2, p2_c0_mesh_glob_temp, iformatted)
    END IF

    !===Start Metis mesh generation=================================================
    ALLOCATE(part(p1_mesh_glob%me))
    WRITE(tit_part,'(i4)') nb_s
    mesh_part_name='mesh_part_S'//trim(adjustl(tit_part))//'.'//trim(adjustl(file_name))
    IF (if_read_partition) THEN
       OPEN(unit=51, file=mesh_part_name, status='unknown', form='formatted')
       READ(51,*) part
       CLOSE(51)
       WRITE(*,*) 'read partition'
    ELSE
       WRITE(*,*) 'create partition'
       CALL part_mesh_m_t_h_phi(nb_s, list_dom_ns, list_dom_temp, list_dom_h, &
            list_dom_phi, p1_mesh_glob, list_inter, part, my_periodic)
       IF (petsc_rank==0) THEN
          OPEN(unit=51, file=mesh_part_name, status='replace', form='formatted')
          WRITE(51,*) part
          CLOSE(51)
       END IF
    END IF

    !===Extract local meshes from global meshes
    IF (if_momentum) THEN
       CALL extract_mesh(comm_one_d(1),nb_s,p1_mesh_glob,part,list_dom_ns,pp_mesh_glob,pp_mesh)
       CALL extract_mesh(comm_one_d(1),nb_s,p2_mesh_glob,part,list_dom_ns,vv_mesh_glob,vv_mesh)

       ALLOCATE(comm_one_d_ns(2))
       comm_one_d_ns(2) = comm_one_d(2) 
       CALL mpi_comm_rank(comm_one_d(1),rank_s,code)
       IF (pp_mesh%me/=0) THEN
          CALL mpi_comm_split(comm_one_d(1),1,rank_s,comm_one_d_ns(1),code)
       ELSE
          CALL mpi_comm_split(comm_one_d(1),mpi_undefined,rank_s,comm_one_d_ns(1),code)
       END IF
    END IF

    IF (if_energy) THEN
       CALL extract_mesh(comm_one_d(1),nb_s,p2_c0_mesh_glob_temp,part,list_dom_temp,temp_mesh_glob,temp_mesh)
       ALLOCATE(comm_one_d_temp(2))
       CALL mpi_comm_dup(comm_one_d(2), comm_one_d_temp(2), code)
       CALL mpi_comm_rank(comm_one_d(1),rank_s,code)
       IF (temp_mesh%me/=0) THEN
          CALL mpi_comm_split(comm_one_d(1),1,rank_s,comm_one_d_temp(1),code)
       ELSE
          CALL mpi_comm_split(comm_one_d(1),mpi_undefined,rank_s,comm_one_d_temp(1),code)
       END IF
    END IF

    IF (if_induction) THEN
       IF (type_fe_h==1) THEN
          CALL extract_mesh(comm_one_d(1),nb_s,p1_mesh_glob,part,list_dom_h,h_mesh_glob,h_mesh)
       ELSE
          CALL extract_mesh(comm_one_d(1),nb_s,p2_mesh_glob,part,list_dom_h,h_mesh_glob,h_mesh)
       END IF
       IF (type_fe_phi==1) THEN
          CALL extract_mesh(comm_one_d(1),nb_s,p1_mesh_glob,part,list_dom_phi,phi_mesh_glob,phi_mesh)
       ELSE
          CALL extract_mesh(comm_one_d(1),nb_s,p2_mesh_glob,part,list_dom_phi,phi_mesh_glob,phi_mesh)
       END IF
    END IF

    !===Cleanup
    CALL free_mesh(p1_mesh_glob)
    CALL free_mesh(p2_mesh_glob)
    IF (if_energy) THEN
       DEALLOCATE(list_inter_temp)
       CALL free_mesh(p2_c0_mesh_glob_temp)
    END IF

    m_max_c = nb_mode/nb_f

    !===Load meshes for monoproc
    IF (if_momentum) THEN
       CALL free_mesh(vv_mesh_glob)
       CALL free_mesh(pp_mesh_glob)
       CALL load_mesh_free_format(directory_m, file_name_m, list_dom_ns,  2, vv_mesh_glob, is_form_m)
       CALL load_mesh_free_format(directory_m, file_name_m, list_dom_ns,  1, pp_mesh_glob, is_form_m)
       IF (check_plt) THEN
          CALL plot_const_p1_label(vv_mesh_glob%jj, vv_mesh_glob%rr, 1.d0*vv_mesh_glob%i_d, 'vv.plt')
       END IF
    END IF
    IF (if_energy) THEN
       CALL free_mesh(temp_mesh_glob)
       CALL load_mesh_free_format(directory_m, file_name_m, list_dom_temp,  2, temp_mesh_glob, is_form_m)
       IF (check_plt) THEN
          CALL plot_const_p1_label(temp_mesh_glob%jj, temp_mesh_glob%rr, 1.d0*temp_mesh_glob%i_d, 'temp.plt')
       END IF
    END IF
    IF (if_induction) THEN
       CALL free_mesh(h_mesh_glob)
       CALL free_mesh(phi_mesh_glob)
       CALL load_dg_mesh_free_format(directory_m, file_name_m, list_dom_h, list_inter_mu,  type_fe_h, h_mesh_glob, is_form_m)
       CALL load_mesh_free_format(directory_m, file_name_m, list_dom_phi, type_fe_phi, phi_mesh_glob, is_form_m)
       IF (check_plt) THEN
          CALL plot_const_p1_label(h_mesh_glob%jj, h_mesh_glob%rr, 1.d0*h_mesh_glob%i_d, 'HH.plt')
          CALL plot_const_p1_label(phi_mesh_glob%jj, phi_mesh_glob%rr, 1.d0*phi_mesh_glob%i_d, 'phi.plt')
       END IF
    END IF

    !===Array allocation
    IF (if_momentum) THEN
       ALLOCATE(un_glob(vv_mesh_glob%np, 6, m_max_c))
       ALLOCATE(un_m1_glob(vv_mesh_glob%np, 6, m_max_c))
       ALLOCATE(pn_glob(pp_mesh_glob%np, 2, m_max_c))
       ALLOCATE(pn_m1_glob(pp_mesh_glob%np, 2, m_max_c))
       ALLOCATE(incpn_m1_glob(pp_mesh_glob%np, 2, m_max_c))
       ALLOCATE(incpn_glob(pp_mesh_glob%np, 2, m_max_c))
       ALLOCATE(un(vv_mesh%np, 6, m_max_c))
       ALLOCATE(un_m1(vv_mesh%np, 6, m_max_c))
       ALLOCATE(pn(pp_mesh%np, 2, m_max_c))
       ALLOCATE(pn_m1(pp_mesh%np, 2, m_max_c))
       ALLOCATE(incpn_m1(pp_mesh%np, 2, m_max_c))
       ALLOCATE(incpn(pp_mesh%np, 2, m_max_c))
       un_glob       = 0.d0
       un_m1_glob    = 0.d0
       pn_glob       = 0.d0
       pn_m1_glob    = 0.d0
       incpn_m1_glob = 0.d0
       incpn_glob    = 0.d0
       un      = 0.d0
       un_m1   = 0.d0
       pn      = 0.d0
       pn_m1   = 0.d0
       incpn_m1= 0.d0
       incpn   = 0.d0
       IF (if_mass) THEN
          ALLOCATE(level_setn_m1_glob(nb_fluid-1,pp_mesh_glob%np, 2, m_max_c))
          ALLOCATE(level_setn_glob(nb_fluid-1,pp_mesh_glob%np, 2, m_max_c))
          ALLOCATE(level_setn_m1(nb_fluid-1,pp_mesh%np, 2, m_max_c))
          ALLOCATE(level_setn(nb_fluid-1,pp_mesh%np, 2, m_max_c))
          level_setn_m1_glob = 0.d0
          level_setn_glob    = 0.d0
          level_setn_m1      = 0.d0
          level_setn         = 0.d0
       END IF
    END IF

    IF (if_energy) THEN
       ALLOCATE(tempn_m1_glob(temp_mesh_glob%np, 2, m_max_c))
       ALLOCATE(tempn_glob(temp_mesh_glob%np, 2, m_max_c))
       ALLOCATE(tempn_m1(temp_mesh%np, 2, m_max_c))
       ALLOCATE(tempn(temp_mesh%np, 2, m_max_c))
       tempn_m1_glob = 0.d0
       tempn_glob    = 0.d0
       tempn_m1      = 0.d0
       tempn         = 0.d0
    END IF

    IF (if_induction) THEN
       ALLOCATE(hn1(h_mesh%np,  6,  m_max_c))
       ALLOCATE(hn(h_mesh%np,  6,  m_max_c))
       ALLOCATE(bn1(h_mesh%np,  6,  m_max_c))
       ALLOCATE(bn(h_mesh%np,  6,  m_max_c))
       ALLOCATE(phin1(phi_mesh%np,2,  m_max_c))
       ALLOCATE(phin(phi_mesh%np,2,  m_max_c))
       ALLOCATE(hn1_glob(h_mesh_glob%np,  6,  m_max_c))
       ALLOCATE(hn_glob(h_mesh_glob%np,  6,  m_max_c))
       ALLOCATE(bn1_glob(h_mesh_glob%np,  6,  m_max_c))
       ALLOCATE(bn_glob(h_mesh_glob%np,  6,  m_max_c))
       ALLOCATE(phin1_glob(phi_mesh_glob%np,2,  m_max_c))
       ALLOCATE(phin_glob(phi_mesh_glob%np,2,  m_max_c))
       hn1        = 0.d0
       hn         = 0.d0
       bn1        = 0.d0
       bn         = 0.d0
       phin1      = 0.d0
       phin       = 0.d0
       hn1_glob   = 0.d0
       hn_glob    = 0.d0
       bn1_glob   = 0.d0
       bn_glob    = 0.d0
       phin1_glob = 0.d0
       phin_glob  = 0.d0
    END IF

    !===Pointers
    IF (is_in) THEN
       mono_in = .false.
       mono_out = .true.
       IF (if_induction) THEN
          hn_in     => hn
          hn1_in    => hn1
          bn_in     => bn
          bn1_in    => bn1
          phin_in   => phin
          phin1_in  => phin1
          h_mesh_in => h_mesh
          phi_mesh_in => phi_mesh
          hn_out     => hn_glob
          hn1_out    => hn1_glob
          bn_out     => bn_glob
          bn1_out    => bn1_glob
          phin_out   => phin_glob
          phin1_out  => phin1_glob
          h_mesh_out => h_mesh_glob
          phi_mesh_out => phi_mesh_glob
       END IF
       IF (if_momentum) THEN
          un_in      => un
          un_m1_in   => un_m1
          pn_in      => pn
          pn_m1_in   => pn_m1
          incpn_in   => incpn
          incpn_m1_in=> incpn_m1
          vv_mesh_in => vv_mesh
          pp_mesh_in => pp_mesh
          un_out      => un_glob
          un_m1_out   => un_m1_glob
          pn_out      => pn_glob
          pn_m1_out   => pn_m1_glob
          incpn_out   => incpn_glob
          incpn_m1_out=> incpn_m1_glob
          vv_mesh_out => vv_mesh_glob
          pp_mesh_out => pp_mesh_glob
          IF (if_mass) THEN
             level_setn_in    => level_setn
             level_setn_m1_in => level_setn_m1
             level_setn_out   => level_setn_glob
             level_setn_m1_out=> level_setn_m1_glob
          END IF
       END IF
       IF (if_energy) THEN
          tempn_in    => tempn
          tempn_m1_in => tempn_m1
          tempn_out   => tempn_glob
          tempn_m1_out => tempn_m1_glob
          temp_mesh_in => temp_mesh
          temp_mesh_out => temp_mesh_glob
       END IF

    ELSE
       mono_in = .true.
       mono_out = .false.
       IF (if_induction) THEN
          hn_in     => hn_glob
          hn1_in    => hn1_glob
          bn_in     => bn_glob
          bn1_in    => bn1_glob
          phin_in   => phin_glob
          phin1_in  => phin1_glob
          h_mesh_in => h_mesh_glob
          phi_mesh_in => phi_mesh_glob
          hn_out     => hn
          hn1_out    => hn1
          bn_out     => bn
          bn1_out    => bn1
          phin_out   => phin
          phin1_out  => phin1
          h_mesh_out => h_mesh
          phi_mesh_out => phi_mesh
       END IF
       IF (if_momentum) THEN
          un_in      => un_glob
          un_m1_in   => un_m1_glob
          pn_in      => pn_glob
          pn_m1_in   => pn_m1_glob
          incpn_in   => incpn_glob
          incpn_m1_in=> incpn_m1_glob
          vv_mesh_in => vv_mesh_glob
          pp_mesh_in => pp_mesh_glob
          un_out      => un
          un_m1_out   => un_m1
          pn_out      => pn
          pn_m1_out   => pn_m1
          incpn_out   => incpn
          incpn_m1_out=> incpn_m1
          vv_mesh_out => vv_mesh
          pp_mesh_out => pp_mesh
          IF (if_mass) THEN
             level_setn_in    => level_setn_glob
             level_setn_m1_in => level_setn_m1_glob
             level_setn_out   => level_setn
             level_setn_m1_out=> level_setn_m1
          END IF
       END IF
       IF (if_energy) THEN
          tempn_in    => tempn_glob
          tempn_m1_in => tempn_m1_glob
          tempn_out   => tempn
          tempn_m1_out=> tempn_m1
          temp_mesh_in => temp_mesh_glob
          temp_mesh_out => temp_mesh
       END IF
    END IF

    !===Interpolation for Maxwell
    IF (rw_mxw) THEN
       IF (rank==0) WRITE(*,*) 'Start interpolation Maxwell'
       IF (inter_mesh) THEN
          ALLOCATE(controle_h(h_mesh_out%np), controle_phi(phi_mesh_out%np))
          DO m = index_start, index_start+nb_fic-1
             hn1        = 0.d0
             hn         = 0.d0
             bn1        = 0.d0
             bn         = 0.d0
             phin1      = 0.d0
             phin       = 0.d0
             hn1_glob   = 0.d0
             hn_glob    = 0.d0
             bn1_glob   = 0.d0
             bn_glob    = 0.d0
             phin1_glob = 0.d0
             phin_glob  = 0.d0
             WRITE(tit, '(i3)') m
             lblank = eval_blank(3,tit)
             DO l = 1, lblank - 1
                tit(l:l) = '0'
             END DO

             IF (petsc_rank==0) CALL system('mv suite_maxwell_I'//tit//'.'//old_filename//'suite_maxwell.'//old_filename)
             CALL mpi_barrier( mpi_comm_world, code)
             CALL read_restart_maxwell(comm_one_d, h_mesh_in, phi_mesh_in, time_h, list_mode, hn_in, hn1_in, &
                  bn_in, bn1_in, phin_in, phin1_in, old_filename, interpol=.false., opt_mono=mono_in)

             !===Controle_H and controle_phi are initialized to zero in interp_mesh
             CALL interp_mesh(h_mesh_in, h_mesh_out, hn_in, hn_out, controle_h, type_fe_h)
             CALL interp_mesh(h_mesh_in, h_mesh_out, hn1_in, hn1_out, controle_h, type_fe_h)
             CALL interp_mesh(h_mesh_in, h_mesh_out, bn_in, bn_out, controle_h, type_fe_h)
             CALL interp_mesh(h_mesh_in, h_mesh_out, bn1_in, bn1_out, controle_h, type_fe_h)
             CALL interp_mesh(phi_mesh_in, phi_mesh_out, phin_in, phin_out, controle_phi, type_fe_phi)
             CALL interp_mesh(phi_mesh_in, phi_mesh_out, phin1_in, phin1_out, controle_phi, type_fe_phi)

             IF ( min(minval(controle_h),minval(controle_phi)) == 0) THEN
                CALL error_petsc('certains points non trouve H/phi 2')
             END IF

             CALL write_restart_maxwell(comm_one_d, h_mesh_out, phi_mesh_out, time_h, list_mode, hn_out, hn1_out, &
                  bn_out, bn1_out, phin_out, phin1_out, new_filename, m, 1, opt_mono=mono_out)
             CALL mpi_comm_rank(comm_one_d(1), rank_s, ierr)
          END DO
          DEALLOCATE(controle_h, controle_phi)
       ELSE
          ALLOCATE(l_t_g_h(h_mesh%np), l_t_g_phi(phi_mesh%np))
          l_t_g_h = 0
          l_t_g_phi = 0
          CALL loc_to_glob(h_mesh, h_mesh_glob, l_t_g_h)
          CALL loc_to_glob(phi_mesh, phi_mesh_glob, l_t_g_phi)
          DO m = index_start, index_start+nb_fic-1
             hn1        = 0.d0
             hn         = 0.d0
             bn1        = 0.d0
             bn         = 0.d0
             phin1      = 0.d0
             phin       = 0.d0
             hn1_glob   = 0.d0
             hn_glob    = 0.d0
             bn1_glob   = 0.d0
             bn_glob    = 0.d0
             phin1_glob = 0.d0
             phin_glob  = 0.d0
             WRITE(tit, '(i3)') m
             lblank = eval_blank(3,tit)
             DO l = 1, lblank - 1
                tit(l:l) = '0'
             END DO

             IF (is_in) THEN
                CALL mpi_comm_rank(comm_one_d(1), rank_s, ierr)
                WRITE(tit_s,'(i3)') rank_s
                lblank = eval_blank(3,tit_s)
                DO l = 1, lblank - 1
                   tit_s(l:l) = '0'
                END DO
                CALL system('mv suite_maxwell_S'//tit_s//'_I'//tit//'.'//old_filename//'suite_maxwell_S'//tit_s//'.'//old_filename)
             ELSE
                IF (petsc_rank==0) CALL system('mv suite_maxwell_I'//tit//'.'//old_filename//'suite_maxwell.'//old_filename)
                CALL mpi_barrier( mpi_comm_world, code)
             END IF
             CALL read_restart_maxwell(comm_one_d, h_mesh_in, phi_mesh_in, time_h, list_mode, hn_in, hn1_in, &
                  bn_in, bn1_in, phin_in, phin1_in, old_filename, interpol=.false., opt_mono=mono_in)

             CALL inter_mesh_loc_to_glob(h_mesh_in, h_mesh_out, hn_in, hn_out, l_t_g_h, is_in, comm_one_d(1))
             CALL inter_mesh_loc_to_glob(h_mesh_in, h_mesh_out, hn1_in, hn1_out, l_t_g_h, is_in, comm_one_d(1))
             CALL inter_mesh_loc_to_glob(h_mesh_in, h_mesh_out, bn_in, bn_out, l_t_g_h, is_in, comm_one_d(1))
             CALL inter_mesh_loc_to_glob(h_mesh_in, h_mesh_out, bn1_in, bn1_out, l_t_g_h, is_in, comm_one_d(1))
             CALL inter_mesh_loc_to_glob(phi_mesh_in, phi_mesh_out, phin_in, phin_out, l_t_g_phi, is_in, comm_one_d(1))
             CALL inter_mesh_loc_to_glob(phi_mesh_in, phi_mesh_out, phin1_in, phin1_out, l_t_g_phi, is_in, comm_one_d(1))

             CALL write_restart_maxwell(comm_one_d, h_mesh_out, phi_mesh_out, time_h, list_mode, hn_out, hn1_out, &
                  bn_out, bn1_out, phin_out, phin1_out, new_filename, m, 1, opt_mono=mono_out)
             CALL mpi_comm_rank(comm_one_d(1), rank_s, ierr)
          END DO

       END IF

       IF (check_plt) THEN
          IF (h_mesh%me /= 0) THEN
             CALL mpi_comm_rank(comm_one_d(1), rang_s, ierr)
             WRITE(tit_s,'(i3)') rang_s
             lblank = eval_blank(3,tit_s)
             DO l = 1, lblank - 1
                tit_s(l:l) = '0'
             END DO
             DO i = 1, SIZE(list_mode)
                WRITE(tit_m,'(i3)') list_mode(i)
                lblank = eval_blank(3,tit_m)
                DO l = 1, lblank - 1
                   tit_m(l:l) = '0'
                END DO
                CALL plot_scalar_field(h_mesh%jj, h_mesh%rr, hn(:,1,i), 'H_r_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
                CALL plot_scalar_field(h_mesh%jj, h_mesh%rr, hn(:,2,i), 'H_r_sin_m='//tit_m//'_'//tit_s//'_999.plt' )
                CALL plot_scalar_field(h_mesh%jj, h_mesh%rr, hn(:,3,i), 'H_t_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
                CALL plot_scalar_field(h_mesh%jj, h_mesh%rr, hn(:,4,i), 'H_t_sin_m='//tit_m//'_'//tit_s//'_999.plt' )
                CALL plot_scalar_field(h_mesh%jj, h_mesh%rr, hn(:,5,i), 'H_z_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
                CALL plot_scalar_field(h_mesh%jj, h_mesh%rr, hn(:,6,i), 'H_z_sin_m='//tit_m//'_'//tit_s//'_999.plt' )

             END DO

             IF (rang_s == 0) THEN
                DO i = 1, SIZE(list_mode)
                   WRITE(tit_m,'(i3)') list_mode(i)
                   lblank = eval_blank(3,tit_m)
                   DO l = 1, lblank - 1
                      tit_m(l:l) = '0'
                   END DO
                   CALL plot_scalar_field(h_mesh_glob%jj, h_mesh_glob%rr, hn_glob(:,1,i), 'gH_r_cos_m='//tit_m//'_999.plt' )
                   CALL plot_scalar_field(h_mesh_glob%jj, h_mesh_glob%rr, hn_glob(:,2,i), 'gH_r_sin_m='//tit_m//'_999.plt' )
                   CALL plot_scalar_field(h_mesh_glob%jj, h_mesh_glob%rr, hn_glob(:,3,i), 'gH_t_cos_m='//tit_m//'_999.plt' )
                   CALL plot_scalar_field(h_mesh_glob%jj, h_mesh_glob%rr, hn_glob(:,4,i), 'gH_t_sin_m='//tit_m//'_999.plt' )
                   CALL plot_scalar_field(h_mesh_glob%jj, h_mesh_glob%rr, hn_glob(:,5,i), 'gH_z_cos_m='//tit_m//'_999.plt' )     
                   CALL plot_scalar_field(h_mesh_glob%jj, h_mesh_glob%rr, hn_glob(:,6,i), 'gH_z_sin_m='//tit_m//'_999.plt' )     
                END DO
             END IF
          END IF
       END IF
       IF (rank==0) WRITE(*,*) 'End interpolation Maxwell'

    END IF

    !===Interpolation for NS
    IF (rw_ns) THEN
       IF (rank==0) WRITE(*,*) 'Start interpolation NS'

       IF (inter_mesh) THEN

          ALLOCATE(controle_vv(vv_mesh_out%np), controle_pp(pp_mesh_out%np))
          DO m = index_start, index_start+nb_fic-1
             un_glob       = 0.d0
             un_m1_glob    = 0.d0
             pn_glob       = 0.d0
             pn_m1_glob    = 0.d0
             incpn_m1_glob = 0.d0
             incpn_glob    = 0.d0
             un      = 0.d0
             un_m1   = 0.d0
             pn      = 0.d0
             pn_m1   = 0.d0
             incpn_m1= 0.d0
             incpn   = 0.d0
             IF (if_level_set) THEN
                level_setn_m1_glob = 0.d0
                level_setn_glob    = 0.d0
                level_setn_m1      = 0.d0
                level_setn         = 0.d0
             END IF

             WRITE(tit, '(i3)') m
             lblank = eval_blank(3,tit)
             DO l = 1, lblank - 1
                tit(l:l) = '0'
             END DO

             IF (petsc_rank==0) THEN
                CALL system('mv suite_ns_I'//tit//'.'//old_filename//'suite_ns.'//old_filename)
             END IF
             CALL mpi_barrier( mpi_comm_world, code)

             IF (vv_mesh%me/=0) THEN
                IF (if_level_set)THEN 
                   CALL read_restart_ns(comm_one_d_ns, time_u, list_mode, un_in, un_m1_in, &
                        pn_in, pn_m1_in, incpn_in, incpn_m1_in, old_filename,opt_level_set=level_setn_in, &
                        opt_level_set_m1=level_setn_m1_in, opt_max_vel=max_vel, opt_mono = mono_in)

                ELSE
                   CALL read_restart_ns(comm_one_d_ns, time_u, list_mode, un_in, un_m1_in, &
                        pn_in, pn_m1_in, incpn_in, incpn_m1_in, old_filename, opt_mono = mono_in)
                END IF
             END IF
             controle_vv = 0
             controle_pp = 0
             CALL interp_mesh(vv_mesh_in, vv_mesh_out, un_in, un_out, controle_vv, 2)
             CALL interp_mesh(vv_mesh_in, vv_mesh_out, un_m1_in, un_m1_out, controle_vv, 2)
             CALL interp_mesh(pp_mesh_in, pp_mesh_out, pn_in, pn_out, controle_pp, 1)
             CALL interp_mesh(pp_mesh_in, pp_mesh_out, pn_m1_in, pn_m1_out, controle_pp, 1)
             CALL interp_mesh(pp_mesh_in, pp_mesh_out, incpn_in, incpn_out, controle_pp, 1)
             CALL interp_mesh(pp_mesh_in, pp_mesh_out, incpn_m1_in, incpn_m1_out, controle_pp, 1)
             IF (if_level_set) THEN
                DO k = 1, nb_fluid-1
                   CALL interp_mesh(pp_mesh_in, pp_mesh_out, level_setn_in(k,:,:,:), level_setn_out(k,:,:,:), controle_pp, 1)
                   CALL interp_mesh(pp_mesh_in, pp_mesh_out, level_setn_m1_in(k,:,:,:), level_setn_m1_out(k,:,:,:), controle_pp, 1)
                END DO
             END IF

             IF ( min(minval(controle_vv),minval(controle_pp)) == 0) THEN
                CALL error_petsc('certains points non trouve')
             END IF

             IF (pp_mesh%me /=0) THEN
                IF (if_level_set) THEN
                   CALL write_restart_ns(comm_one_d_ns, vv_mesh_out, pp_mesh_out, time_u, list_mode, un_out, un_m1_out, &
                        pn_out, pn_m1_out, incpn_out, incpn_m1_out, new_filename, m, 1, opt_level_set=level_setn_out, &
                        opt_level_set_m1=level_setn_m1_out, opt_max_vel=max_vel, opt_mono=mono_out)
                ELSE
                   CALL write_restart_ns(comm_one_d_ns, vv_mesh_out, pp_mesh_out, time_u, list_mode, un_out, un_m1_out, &
                        pn_out, pn_m1_out, incpn_out, incpn_m1_out, new_filename, m, 1, opt_mono=mono_out)
                END IF
                CALL mpi_comm_rank(comm_one_d_ns(1), rang_ns_s, ierr)
             END IF
          END DO
          DEALLOCATE(controle_vv, controle_pp)

       ELSE

          ALLOCATE(l_t_g_vv(vv_mesh%np), l_t_g_pp(pp_mesh%np))
          l_t_g_vv = 0
          l_t_g_pp = 0
          CALL loc_to_glob(vv_mesh, vv_mesh_glob, l_t_g_vv)
          CALL loc_to_glob(pp_mesh, pp_mesh_glob, l_t_g_pp)
          IF (pp_mesh%me /=0) THEN
             DO m = index_start, index_start+nb_fic-1
                un_glob       = 0.d0
                un_m1_glob    = 0.d0
                pn_glob       = 0.d0
                pn_m1_glob    = 0.d0
                incpn_m1_glob = 0.d0
                incpn_glob    = 0.d0
                un      = 0.d0
                un_m1   = 0.d0
                pn      = 0.d0
                pn_m1   = 0.d0
                incpn_m1= 0.d0
                incpn   = 0.d0
                IF (if_level_set) THEN
                   level_setn_m1_glob = 0.d0
                   level_setn_glob    = 0.d0
                   level_setn_m1      = 0.d0
                   level_setn         = 0.d0
                END IF

                WRITE(tit, '(i3)') m
                lblank = eval_blank(3,tit)
                DO l = 1, lblank - 1
                   tit(l:l) = '0'
                END DO
                IF (is_in) THEN
                   CALL mpi_comm_rank(comm_one_d_ns(1), rang_ns_s, ierr)
                   WRITE(tit_s,'(i3)') rang_ns_s
                   lblank = eval_blank(3,tit_s)
                   DO l = 1, lblank - 1
                      tit_s(l:l) = '0'
                   END DO

                   CALL system('mv suite_ns_S'//tit_s//'_I'//tit//'.'//old_filename//'suite_ns_S'//tit_s//'.'//old_filename)

                ELSE
                   IF (petsc_rank==0) CALL system('mv suite_ns_I'//tit//'.'//old_filename//'suite_ns.'//old_filename)
                   CALL mpi_barrier( mpi_comm_world, code)
                END IF

                IF (if_level_set) THEN
                   CALL read_restart_ns(comm_one_d_ns, time_u, list_mode, un_in, un_m1_in, &
                        pn_in, pn_m1_in, incpn_in, incpn_m1_in, old_filename,  opt_level_set=level_setn_in, &
                        opt_level_set_m1=level_setn_m1_in, opt_max_vel=max_vel, opt_mono = mono_in)
                ELSE
                   CALL read_restart_ns(comm_one_d_ns, time_u, list_mode, un_in, un_m1_in, &
                        pn_in, pn_m1_in, incpn_in, incpn_m1_in, old_filename, opt_mono = mono_in)
                END IF

                CALL inter_mesh_loc_to_glob(vv_mesh_in, vv_mesh_out, un_in, un_out, l_t_g_vv, is_in, &
                     comm_one_d_ns(1))
                CALL inter_mesh_loc_to_glob(vv_mesh_in, vv_mesh_out, un_m1_in, un_m1_out, l_t_g_vv, is_in, &
                     comm_one_d_ns(1))
                CALL inter_mesh_loc_to_glob(pp_mesh_in, pp_mesh_out, pn_in, pn_out, l_t_g_pp, is_in, & 
                     comm_one_d_ns(1))
                CALL inter_mesh_loc_to_glob(pp_mesh_in, pp_mesh_out, pn_m1_in, pn_m1_out, l_t_g_pp, is_in, &
                     comm_one_d_ns(1))
                CALL inter_mesh_loc_to_glob(pp_mesh_in, pp_mesh_out, incpn_in, incpn_out, l_t_g_pp, is_in, &
                     comm_one_d_ns(1))
                CALL inter_mesh_loc_to_glob(pp_mesh_in, pp_mesh_out, incpn_m1_in, incpn_m1_out, l_t_g_pp, is_in, &
                     comm_one_d_ns(1))
                IF (if_level_set) THEN
                   DO k = 1, nb_fluid-1
                      CALL inter_mesh_loc_to_glob(pp_mesh_in, pp_mesh_out, level_setn_in(k,:,:,:), level_setn_out(k,:,:,:), &
                           l_t_g_pp, is_in, comm_one_d_ns(1))
                      CALL inter_mesh_loc_to_glob(pp_mesh_in, pp_mesh_out, level_setn_m1_in(k,:,:,:), &
                           level_setn_m1_out(k,:,:,:), l_t_g_pp, is_in, comm_one_d_ns(1))
                   END DO
                END IF

                IF (if_level_set) THEN
                   CALL write_restart_ns(comm_one_d_ns, vv_mesh_out, pp_mesh_out, time_u, list_mode, un_out, un_m1_out, &
                        pn_out, pn_m1_out, incpn_out, incpn_m1_out, new_filename, m, 1, opt_level_set=level_setn_out, &
                        opt_level_set_m1=level_setn_m1_out, opt_max_vel=max_vel, opt_mono=mono_out)
                ELSE
                   CALL write_restart_ns(comm_one_d_ns, vv_mesh_out, pp_mesh_out, time_u, list_mode, un_out, un_m1_out, &
                        pn_out, pn_m1_out, incpn_out, incpn_m1_out, new_filename, m, 1, opt_mono=mono_out)
                END IF
                CALL mpi_comm_rank(comm_one_d_ns(1), rang_ns_s, ierr)
             END DO
          END IF
       END IF

       IF (check_plt) THEN
          IF (pp_mesh%me /= 0) THEN
             CALL mpi_comm_rank(comm_one_d_ns(1), rang_ns_s, ierr)
             WRITE(tit_s,'(i3)') rang_ns_s
             lblank = eval_blank(3,tit_s)
             DO l = 1, lblank - 1
                tit_s(l:l) = '0'
             END DO
             DO i = 1, SIZE(list_mode)
                WRITE(tit_m,'(i3)') list_mode(i)
                lblank = eval_blank(3,tit_m)
                DO l = 1, lblank - 1
                   tit_m(l:l) = '0'
                END DO
                CALL plot_scalar_field(vv_mesh%jj, vv_mesh%rr, un(:,1,i), 'u_r_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
                CALL plot_scalar_field(vv_mesh%jj, vv_mesh%rr, un(:,3,i), 'u_t_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
                CALL plot_scalar_field(vv_mesh%jj, vv_mesh%rr, un(:,4,i), 'u_t_sin_m='//tit_m//'_'//tit_s//'_999.plt' )
                CALL plot_scalar_field(vv_mesh%jj, vv_mesh%rr, un(:,5,i), 'u_z_cos_m='//tit_m//'_'//tit_s//'_999.plt' )

             END DO

             IF (rang_ns_s == 0) THEN
                DO i = 1, SIZE(list_mode)
                   WRITE(tit_m,'(i3)') list_mode(i)
                   lblank = eval_blank(3,tit_m)
                   DO l = 1, lblank - 1
                      tit_m(l:l) = '0'
                   END DO
                   CALL plot_scalar_field(vv_mesh_glob%jj, vv_mesh_glob%rr, un_glob(:,1,i), 'gu_r_cos_m='//tit_m//'_999.plt' )
                   CALL plot_scalar_field(vv_mesh_glob%jj, vv_mesh_glob%rr, un_glob(:,3,i), 'gu_t_cos_m='//tit_m//'_999.plt' )
                   CALL plot_scalar_field(vv_mesh_glob%jj, vv_mesh_glob%rr, un_glob(:,4,i), 'gu_t_sin_m='//tit_m//'_999.plt' )
                   CALL plot_scalar_field(vv_mesh_glob%jj, vv_mesh_glob%rr, un_glob(:,5,i), 'gu_z_cos_m='//tit_m//'_999.plt' )     
                END DO
             END IF
          END IF
       END IF
       IF (rank==0) WRITE(*,*) 'End interpolation NS'
    END IF

    !===Interpolation for temperature
    IF (rw_temp) THEN
       IF (rank==0) WRITE(*,*) 'Start interpolation temperature'

       IF (inter_mesh) THEN
          ALLOCATE(controle_temp(temp_mesh_out%np))
          DO m = index_start, index_start+nb_fic-1
             tempn_m1_glob = 0.d0
             tempn_glob    = 0.d0
             tempn_m1      = 0.d0
             tempn         = 0.d0
             WRITE(tit, '(i3)') m
             lblank = eval_blank(3,tit)
             DO l = 1, lblank - 1
                tit(l:l) = '0'
             END DO

             IF (petsc_rank==0) THEN
                CALL system('mv suite_temp_I'//tit//'.'//old_filename//'suite_temp.'//old_filename)
             END IF
             CALL mpi_barrier( mpi_comm_world, code)

             IF (temp_mesh%me/=0) THEN
                CALL read_restart_temp(comm_one_d_temp, time_t, list_mode, tempn_in, tempn_m1_in, old_filename, &
                     opt_mono = mono_in)
             END IF
             controle_temp = 0
             CALL interp_mesh(temp_mesh_in, temp_mesh_out, tempn_in, tempn_out, controle_temp, 2)
             CALL interp_mesh(temp_mesh_in, temp_mesh_out, tempn_m1_in, tempn_m1_out, controle_temp, 2)

             IF (temp_mesh%me /= 0) THEN
                CALL write_restart_temp(comm_one_d_temp, temp_mesh_out, time_t, &
                     list_mode, tempn_out, tempn_m1_out, new_filename, m, 1, opt_mono = mono_out)
                CALL mpi_comm_rank(comm_one_d_temp(1), rang_temp_s, ierr)
             END IF
          END DO
          DEALLOCATE(controle_temp)

       ELSE

          ALLOCATE(l_t_g_temp(temp_mesh%np))
          l_t_g_temp = 0
          CALL loc_to_glob(temp_mesh, temp_mesh_glob, l_t_g_temp)
          IF (temp_mesh%me /=0) THEN
             DO m = index_start, index_start+nb_fic-1
                tempn_m1_glob = 0.d0
                tempn_glob    = 0.d0
                tempn_m1      = 0.d0
                tempn         = 0.d0

                WRITE(tit, '(i3)') m
                lblank = eval_blank(3,tit)
                DO l = 1, lblank - 1
                   tit(l:l) = '0'
                END DO
                IF (is_in) THEN
                   CALL mpi_comm_rank(comm_one_d_temp(1), rang_temp_s, ierr)
                   WRITE(tit_s,'(i3)') rang_temp_s
                   lblank = eval_blank(3,tit_s)
                   DO l = 1, lblank - 1
                      tit_s(l:l) = '0'
                   END DO

                   CALL system('mv suite_temp_S'//tit_s//'_I'//tit//'.'//old_filename//'suite_temp_S'//tit_s//'.'//old_filename)
                ELSE
                   IF (petsc_rank==0) CALL system('mv suite_temp_I'//tit//'.'//old_filename//'suite_temp.'//old_filename)
                   CALL mpi_barrier( mpi_comm_world, code)
                END IF

                CALL read_restart_temp(comm_one_d_temp, time_t, list_mode, tempn_in, tempn_m1_in, old_filename, &
                     opt_mono = mono_in)

                CALL inter_mesh_loc_to_glob(temp_mesh_in, temp_mesh_out, tempn_in, tempn_out, l_t_g_temp, is_in, &
                     comm_one_d_temp(1))
                CALL inter_mesh_loc_to_glob(temp_mesh_in, temp_mesh_out, tempn_m1_in, tempn_m1_out, l_t_g_temp, is_in, &
                     comm_one_d_temp(1))

                CALL write_restart_temp(comm_one_d_temp, temp_mesh_out, time_t, list_mode, tempn_out, tempn_m1_out, &
                     new_filename, m, 1, opt_mono = mono_out)
                CALL mpi_comm_rank(comm_one_d_temp(1), rang_temp_s, ierr)
             END DO
          END IF
       END IF

       IF (check_plt) THEN
          IF (temp_mesh%me /= 0) THEN
             CALL mpi_comm_rank(comm_one_d_temp(1), rang_temp_s, ierr)
             WRITE(tit_s,'(i3)') rang_temp_s
             lblank = eval_blank(3,tit_s)
             DO l = 1, lblank - 1
                tit_s(l:l) = '0'
             END DO
             DO i = 1, SIZE(list_mode)
                WRITE(tit_m,'(i3)') list_mode(i)
                lblank = eval_blank(3,tit_m)
                DO l = 1, lblank - 1
                   tit_m(l:l) = '0'
                END DO
                CALL plot_scalar_field(temp_mesh%jj, temp_mesh%rr, tempn(:,1,i), 'temp_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
                CALL plot_scalar_field(temp_mesh%jj, temp_mesh%rr, tempn(:,2,i), 'temp_sin_m='//tit_m//'_'//tit_s//'_999.plt' )
             END DO
             IF (rang_temp_s == 0) THEN
                DO i = 1, SIZE(list_mode)
                   WRITE(tit_m,'(i3)') list_mode(i)
                   lblank = eval_blank(3,tit_m)
                   DO l = 1, lblank - 1
                      tit_m(l:l) = '0'
                   END DO
                   CALL plot_scalar_field(temp_mesh_glob%jj, temp_mesh_glob%rr, tempn(:,1,i), &
                        'gtemp_cos_m='//tit_m//'_'//tit_s//'_999.plt' )
                   CALL plot_scalar_field(temp_mesh_glob%jj, temp_mesh_glob%rr, tempn(:,2,i), &
                        'gtemp_sin_m='//tit_m//'_'//tit_s//'_999.plt' )
                END DO
             END IF
          END IF
       END IF
       IF (rank==0) WRITE(*,*) 'End interpolation temperature'
    END IF

    IF (is_in) THEN
       IF (petsc_rank==0 .AND. rw_mxw) CALL system('rm -rf suite_maxwell_S*')
       IF (petsc_rank==0 .AND. rw_ns) CALL system('rm -rf suite_ns_S*')
       IF (petsc_rank==0 .AND. rw_temp) CALL system('rm -rf suite_temp_S*')
    ELSE
       IF (petsc_rank==0 .AND. rw_mxw) CALL system('rm -rf suite_maxwell.*')
       IF (petsc_rank==0 .AND. rw_ns) CALL system('rm -rf suite_ns.*')
       IF (petsc_rank==0 .AND. rw_temp) CALL system('rm -rf suite_temp.*')
    END IF

    IF (check_plt .AND. petsc_rank==0) THEN
       CALL system('mkdir -p PLT')
       CALL system('mv *.plt PLT/')
    END IF
    CALL mpi_barrier( mpi_comm_world, code)

    CALL petscfinalize(ierr)

  END SUBROUTINE mesh_interpol

  SUBROUTINE which_pb(is_in)
    USE my_util
    IMPLICIT NONE
    LOGICAL, INTENT(OUT) :: is_in
    CHARACTER(len=200)   :: inline
    CALL getarg(1,inline)
    IF (trim(adjustl(inline)) == 'petsc_nonpetsc') THEN
       is_in = .true.
    ELSE IF (trim(adjustl(inline)) == 'nonpetsc_petsc') THEN
       is_in = .false.
    ELSE
       CALL error_petsc('BUG in which_pb')
    END IF
  END SUBROUTINE which_pb

END MODULE mesh_interpolation