 !***********************************************************************
 !
 !    Copyright (c) 2015, Lawrence Livermore National Security, LLC.
 !                        Produced at the Lawrence Livermore National
 !                        Laboratory.
 !                        Written by Nicolas Schunck, schunck1@llnl.gov
 !
 !    LLNL-CODE-470611 All rights reserved.
 !
 !    Copyright 2011, N. Schunck, J. Dobaczewski, J. McDonnell,
 !                    W. Satula, J.A. Sheikh, A. Staszczak,A
 !                    M. Stoitsov, P. Toivanen
 !    Copyright 2009, J. Dobaczewski, W. Satula, B.G. Carlsson, J. Engel,
 !                    P. Olbratowski, P. Powalowski, M. Sadziak,
 !                    J. Sarich, N. Schunck, A. Staszczak, M. Stoitsov,
 !                    M. Zalewski, H. Zdunczuk
 !    Copyright 2004, 2005, J. Dobaczewski, P. Olbratowski
 !    Copyright 1997, 2000, J. Dobaczewski, J. Dudek
 !
 !    This file is part of HFODD.
 !
 !    HFODD is free software: you can redistribute it and/or modify it
 !    under the terms of the GNU General Public License as published by
 !    the Free Software Foundation, either version 3 of the License, or
 !    (at your option) any later version.
 !
 !    HFODD is distributed in the hope that it will be useful, but
 !    WITHOUT ANY WARRANTY; without even the implied warranty of
 !    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 !    GNU General Public License for more details.
 !
 !    You should have received a copy of the GNU General Public License
 !    along with HFODD. If not, see <http://www.gnu.org/licenses/>.
 !
 !    OUR NOTICE AND TERMS AND CONDITIONS OF THE GNU GENERAL PUBLIC
 !    LICENSE
 !
 !    Our Preamble Notice
 !
 !      A. This notice is required to be provided under our contract
 !         with the U.S. Department of Energy (DOE). This work was
 !         produced at the Lawrence Livermore National Laboratory under
 !         Contract No. DE-AC52-07NA27344 with the DOE.
 !      B. Neither the United States Government nor Lawrence Livermore
 !         National Security, LLC nor any of their employees, makes any
 !         warranty, express or implied, or assumes any liability or
 !         responsibility for the accuracy, completeness, or usefulness
 !         of any information, apparatus, product, or process disclosed,
 !         or represents that its use would not infringe privately-owned
 !         rights.
 !      C. Also, reference herein to any specific commercial products,
 !         process, or services by trade name, trademark, manufacturer
 !         or otherwise does not necessarily constitute or imply its
 !         endorsement, recommendation, or favoring by the United States
 !         Government or Lawrence Livermore National Security, LLC. The
 !         views and opinions of authors expressed herein do not
 !         necessarily state or reflect those of the United States
 !         Government or Lawrence Livermore National Security, LLC, and
 !         shall not be used for advertising or product endorsement
 !         purposes.
 !
 !    The precise terms and conditions for copying, distribution and
 !    modification are contained in the file COPYING.
 !
 !***********************************************************************

 !----------------------------------------------------------------------!
 !                                                                      !
 !                MPI Management Environment for HFODD                  !
 !                                                                      !
 ! Authors: N. Schunck, ORNL                                            !
 !                                                                      !
 ! This module is/will be in charge of:                                 !
 !    - defining the list of tasks to be run based on sequential and    !
 !      parallel inputs read in hfodd_mpiio.f90                         !
 !    - managing the load balancing of the tasks                        !
 !                                                                      !
 !                          First included in official release v249e    !
 !                                                                      !
 !----------------------------------------------------------------------!

      SUBROUTINE Set_version_hfodd_mpimanager
      !
      USE hfodd_sizes
      !
      CHARACTER NAMMOD*16
      INTEGER   MODUVE,MODSET
      !
      COMMON /VERMOD/ NAMMOD(NDMODU),MODUVE(NDMODU),MODSET(NDMODU)
      !
      NAMMOD(7) = 'hfodd_mpimanager'
      MODUVE(7) = 6
      MODSET(7) = 7
      !
      END SUBROUTINE Set_version_hfodd_mpimanager
      !
 Module hfodd_mpimanager

   Use hfodd_mpiio

   Implicit None

   ! Public
   Integer :: NOFJOB

   ! Private
   Logical :: filled

   Integer :: debug = 0

   Integer :: numberQ=4,numberQ_restart=4,numberQ_orig=4,       &
              numberConstraints_orig=1,neck_numberPoints_orig=1
   Integer, Allocatable :: IZ_JOB(:),IN_JOB(:),indexQjob(:),indexForce(:),indexTemperature(:), &
                           indexRestartJob(:),indexNeck(:)
   Integer, Allocatable :: indexSizeJob(:),indexStatesJob(:),indexLengthJob(:),idxDeformJob(:),    &
                           basis_size(:),basis_states(:),flex_lambda(:),flex_miu(:),index_file(:), &
                           mapping_qlm(:),mapping_qlm_new(:),mapping_neck(:)
   Integer, Allocatable :: lambdaJob(:,:),miuJob(:,:),lambdaJob_restart(:,:),miuJob_restart(:,:)
   Integer, Allocatable :: table_fichiers(:)

   Real :: neck_step=0.0,neck_step_restart
   Real, Allocatable :: XIZJOB(:),XINJOB(:)
   Real, Allocatable :: flex_qlm(:),basis_freq(:),basis_defo(:), &
                        neck_value(:),neck_value_restart(:)
   Real, Allocatable :: qValueJob(:,:),qValueJob_restart(:,:)

   Character(Len=68), Allocatable :: fichier(:),fichier_tho(:)
   Character(Len=68), Allocatable :: fichier_rec_old(:),fichier_rec_new(:),&
                                     fichier_tho_old(:),fichier_tho_new(:),&
                                     fichier_lic_old(:),fichier_lic_new(:)
   Character(Len=68), Allocatable :: fichier_rec_tmp(:),&
                                     fichier_tho_tmp(:),&
                                     fichier_lic_tmp(:)

 Contains

   !---------------------------------------------------------------------!
   !  Defines default deformation grid for 0<= modeMPI_deformation <=2   !
   !  Outputs: - numberQ, total number of points on the mesh             !
   !           - lambdaJob(:), miuJob(:), qValueJob(:), vectors storing  !
   !             the lambda, miu and q(lambda,miu) values at each point  !
   !             of the mesh                                             !
   !           - numberQ_orig, numberConstraints_orig store the sizes of !
   !             the mesh (total size and number of constraints, resp.)  !
   !             for use when modeMPI_deformation = 2                    !
   !---------------------------------------------------------------------!
   Subroutine define_gridQlm()

     Integer :: iq,nQtotal,nQfactor,ii,lambda,miu,nQ,mpi_rank,mpi_err
     Integer, Allocatable :: nQmodul(:)
     Real, Allocatable :: deltaQ(:)

     Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err)

     ! Determine total number of points on the mesh: number of collective
     ! variables 'numberConstraints' times the number of points for each
     numberQ = 1
     Do iq=1,numberConstraints
        numberQ = numberQ*qStepNumber(iq)
     End Do
     If(debug.Ge.1 .And. mpi_rank.Eq.0) &
        Write(6,'("numberQ = ",i10," numberConstraints=",i10)') numberQ,numberConstraints

     ! Storing characteristics of the default mesh for future use (if
     ! modeMPI_deformation = 2)
     numberQ_orig = numberQ
     numberConstraints_orig = numberConstraints

     ! Getting the meshsize for each collective variable
     Allocate(deltaQ(1:numberConstraints),nQmodul(0:numberConstraints))
     Do iq=1,numberConstraints
        If(qstepNumber(iq) > 1 ) Then
            deltaQ(iq) = (qFinValue(iq) - qinitValue(iq)) / Real(qstepNumber(iq) - 1)
        Else
            deltaQ(iq) = 0.0d0
        End If
     End Do

     ! Defining the mesh: lambdaJob(ii,iq), miuJob(ii,iq) and qValueJob(ii,iq) give,
     ! respectively, the value of lambda, miu and q(lambda,miu) at point ii of the
     ! total mesh and for the collective variable iq
     Allocate(qValueJob(1:numberQ,1:numberConstraints))
     Allocate(lambdaJob(1:numberQ,1:numberConstraints),miuJob(1:numberQ,1:numberConstraints))

     nQtotal=1; nQmodul(0)=1; nQfactor=0
     Do iq=1,numberConstraints

        nQ=qStepNumber(iq)
        lambda=Abs(qLambda(iq)); miu=qMiu(iq)
        nQtotal=nQtotal*nQ; nQmodul(iq)=numberQ/nQtotal

        Do ii=1,numberQ
           lambdaJob(ii,iq)=lambda
           miuJob(ii,iq)=miu
           qValueJob(ii,iq)=qinitValue(iq) &
                           + deltaQ(iq) * ((ii-1)/nQmodul(iq) - nQfactor*((ii-1)/nQmodul(iq-1))*nQ)
        End Do

        nQfactor = 1

     End Do

     Deallocate(deltaQ,nQmodul)

   End Subroutine define_gridQlm

   !---------------------------------------------------------------------!
   !  Defines restart deformation grid for modeMPI_deformation = 2, 3    !
   !  Outputs: - numberQ_restart, total number of points on the restart  !
   !             mesh                                                    !
   !           - lambdaJob_restart(:), miuJob_restart(:), and            !
   !             qValueJob_restart(:), vectors storing the lambda, miu   !
   !             and q(lambda,miu) values at each point of the restart   !
   !             mesh                                                    !
   !---------------------------------------------------------------------!
   Subroutine define_gridQlm_new()

     Integer :: iq,nQtotal,nQfactor,ii,lambda,miu,nQ,mpi_rank,mpi_err
     Integer, Allocatable :: nQmodul(:)
     Real, Allocatable :: deltaQ(:)

     Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err)

     ! Determine total number of points on the restart mesh
     numberQ_restart = 1
     Do iq=1,numberConstraints_restart
        numberQ_restart = numberQ_restart * qStepNumber_restart(iq)
     End Do
     If(debug.Ge.1 .And. mpi_rank.Eq.0) &
        Write(6,'("numberQ_restart = ",i10," numberConstraints_restart=",i10)') &
                   numberQ_restart,numberConstraints_restart

     ! Getting the meshsize for each collective variable on the restart mesh
     Allocate(deltaQ(1:numberConstraints_restart),nQmodul(0:numberConstraints_restart))
     Do iq=1,numberConstraints_restart
        If(qstepNumber_restart(iq)>1) Then
           deltaQ(iq) = (qFinValue_restart(iq) - qinitValue_restart(iq)) / Real(qstepNumber_restart(iq) - 1)
        Else
           deltaQ(iq) = 0.0D0
        End If
     End Do

     ! Defining the restart mesh, see subroutine define_gridQlm for explanations
     Allocate(qValueJob_restart(1:numberQ_restart,1:numberConstraints_restart))
     Allocate(lambdaJob_restart(1:numberQ_restart,1:numberConstraints_restart),&
              miuJob_restart(1:numberQ_restart,1:numberConstraints_restart))

     nQtotal=1; nQmodul(0)=1; nQfactor=0
     Do iq=1,numberConstraints_restart

        nQ=qStepNumber_restart(iq)
        lambda=Abs(qLambda_restart(iq)); miu=qMiu_restart(iq)
        nQtotal=nQtotal*nQ; nQmodul(iq)=numberQ_restart / nQtotal

        Do ii=1,numberQ_restart
           lambdaJob_restart(ii,iq)=lambda
           miuJob_restart(ii,iq)=miu
           qValueJob_restart(ii,iq)=qinitValue_restart(iq) &
                                   + deltaQ(iq) * ( (ii-1)/nQmodul(iq) - nQfactor*((ii-1)/nQmodul(iq-1))*nQ )
        End Do

        nQfactor = 1

     End Do

     Deallocate(deltaQ,nQmodul)

   End Subroutine define_gridQlm_new

   !---------------------------------------------------------------------!
   !  Defines the deformation path for modeMPI_deformation = 4. The full !
   !  deformation path is stored in the file hfodd_path.d. It is read in !
   !  this routine by processor 0 and broadcast to all processors.       !
   !  Outputs: - numberQ_restart, total number of points in the path     !
   !           - lambdaJob_restart(:), miuJob_restart(:), and            !
   !             qValueJob_restart(:), vectors storing the lambda, miu   !
   !             and q(lambda,miu) values at each point along the path   !
   !---------------------------------------------------------------------!
   Subroutine define_pathQlm_new()

     Integer, Parameter :: ndatin=84
     Integer :: i,lines,iline,number_cons,ierr,mpi_rank,mpi_err
     Character(Len=80) :: filein

     ! Getting the rank of the current process
     Call mpi_comm_rank(MPI_COMM_WORLD, mpi_rank, mpi_err)

     If(mpi_rank.Eq.0) Then

        ! open the file
        filein = "hfodd_path.d"
        Open(ndatin,file=filein,status='old',form='formatted',iostat=ierr)
        If(ierr.Ne.0) Then
           Write(6,'("Error in opening file : ",a80)') filein
           Stop 'Error in define_pathQlm_new - I/O'
        End if

        ! Read the number of constraints used to define a trajectory
        Read(ndatin,*) number_cons, lines
        If(debug.Ge.1) Write(6,'("Number of constraints: ",i10," Number of lines: ",i10)') number_cons,lines
        numberConstraints_restart=number_cons
        numberQ_restart=lines
        ! Read trajectory in N-d space (N=number_cons)
        Allocate(qValueJob_restart(1:numberQ_restart,1:numberConstraints_restart))
        Allocate(lambdaJob_restart(1:numberQ_restart,1:numberConstraints_restart),&
                 miuJob_restart(1:numberQ_restart,1:numberConstraints_restart))
        iline=0
        Do While(iline < lines)
           iline=iline+1
           Read(ndatin,*) (lambdaJob_restart(iline,i),miuJob_restart(iline,i),qValueJob_restart(iline,i),&
                           i=1,number_cons)
        End Do
        ! Close file
        Close(ndatin)

     End If

     ! Broadcast number of constraints and number of lines
     Call MPI_Bcast(numberConstraints_restart,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err)
     Call MPI_Bcast(numberQ_restart,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err)
     ! Allocate space for vectors
     If(mpi_rank.Ne.0) Then
        Allocate(qValueJob_restart(1:numberQ_restart,1:numberConstraints_restart))
        Allocate(lambdaJob_restart(1:numberQ_restart,1:numberConstraints_restart),&
                 miuJob_restart(1:numberQ_restart,1:numberConstraints_restart))
     End If

     ! Broadcast lambda values                                  ,
     Call MPI_Bcast(lambdaJob_restart,numberConstraints_restart*numberQ_restart,MPI_INTEGER,0,&
                                                                                MPI_COMM_WORLD,mpi_err)
     ! Broadcast miu values                                     ,
     Call MPI_Bcast(miuJob_restart,numberConstraints_restart*numberQ_restart,MPI_INTEGER,0,&
                                                                             MPI_COMM_WORLD,mpi_err)
     ! Broadcast qlm values
     Call MPI_Bcast(qValueJob_restart,numberConstraints_restart*numberQ_restart,MPI_DOUBLE_PRECISION,0,&
                                                                                MPI_COMM_WORLD,mpi_err)
     If(debug.Ge.1) Write(6,'("In define_pathQLM_new - mpi_rank = ",i10, &
                            & " has received numberConstraints_restart = ",i10)') &
                          mpi_rank, numberConstraints_restart

   End Subroutine define_pathQlm_new

   !---------------------------------------------------------------------!
   !  Defines default grid for the constraint on the size of the neck QN !
   !  in the case modeMPI_neck >= 1                                      !
   !  Outputs: - neck_value(:), vector storing the values of QN          !
   !           - neck_numberPoints_orig stores the number of QN values   !
   !             for use when modeMPI_neck = 2                           !
   !---------------------------------------------------------------------!
   Subroutine define_gridNeck()

     Integer :: i_neck

     neck_numberPoints_orig = neck_numberPoints

     Allocate(neck_value(1:neck_numberPoints))
     If(neck_numberPoints>1) Then
         neck_step = (neck_finValue - neck_initValue) / Real(neck_numberPoints - 1)
         do i_neck=1,neck_numberPoints
            neck_value(i_neck) = neck_initValue + Real(i_neck - 1)*neck_step
         End do
     else
         neck_value(1) = neck_initValue
     End if

   End Subroutine define_gridNeck

   !---------------------------------------------------------------------!
   !  Defines restart grid for the constraint on the size of the neck QN !
   !  in the case modeMPI_neck = 2                                       !
   !  Outputs: - neck_value_restart(:), vector storing the values of QN  !
   !             on the restart grid                                     !
   !---------------------------------------------------------------------!
   Subroutine define_gridNeck_new()

     Integer :: i_neck

     Allocate(neck_value_restart(1:neck_numberPoints_restart))
     If(neck_numberPoints_restart>1) Then
         neck_step_restart = (neck_finValue_restart - neck_initValue_restart) &
                          / Real(neck_numberPoints_restart - 1)
         do i_neck=1,neck_numberPoints_restart
            neck_value_restart(i_neck) = neck_initValue_restart + Real(i_neck - 1)*neck_step_restart
         End do
     else
         neck_value_restart(1) = neck_initValue_restart
     End if

   End Subroutine define_gridNeck_new

   !---------------------------------------------------------------------!
   !  Defines default deformation grid for modeMPI_deformation >= 3. In  !
   !  this mode, the initial mesh is extracted on-the-fly from a list of !
   !  HFODD record files. As a result, the mesh can be irregular and     !
   !  have different number of constraints from the restart mesh. The    !
   !  reading of the HFODD record files is done in parallel to handle    !
   !  large file sets (and reducing the overhead of reading HFODD binary !
   !  files).                                                            !
   !  ATTENTION: We make the assumption that all files in the list have  !
   !             the same number (=numberConstraints) and type (qLambda, !
   !             qMiu arrays) of constraints. This information is extrac-!
   !             ted from the MULTICONS keyword. It is necessary to spe- !
   !             cify what the constraints are, even though the actual   !
   !             mesh will be discarded for these two modes.             !
   !  Outputs: - numberQ, total number of points on the mesh             !
   !           - lambdaJob(:), miuJob(:), qValueJob(:), vectors storing  !
   !             the lambda, miu and q(lambda,miu) values at each point  !
   !             of the mesh                                             !
   !           - numberQ_orig, numberConstraints_orig store the sizes of !
   !             the mesh (total size and number of constraints, resp.)  !
   !           - neck_numberPoints_orig is set to 1, as this mode is not !
   !             yet compatible with any constraint on the neck.         !
   !---------------------------------------------------------------------!
   Subroutine define_restartPointsQlm(NUMITE,NXHERM,NYHERM,NZHERM,IPCONT,ILCONT,IACONT, &
                                      IMCONT,IRCONT,IRENMA,IRENIN,LIPKIN,LIPKIP,REFERN, &
                                      REFERP,REDELN,REDELP,REFE2N,REFE2P,IFIBLN,INIBLN, &
                                      IFIBLP,INIBLP,NMUMAX,NSIMAX,NMUCON,ISHIFT,IF_THO, &
                                                                                IERROR)

     Integer, INTENT(IN) :: NUMITE,NXHERM,NYHERM,NZHERM,IPCONT,ILCONT,IACONT,IMCONT, &
                            IRCONT,IRENMA,IRENIN,LIPKIN,LIPKIP,IFIBLN,INIBLN,IFIBLP, &
                            INIBLP,NMUMAX,NSIMAX,NMUCON,ISHIFT,IF_THO
     Integer, INTENT(INOUT) :: IERROR
     Real, INTENT(INOUT) :: REFERN,REFERP,REDELN,REDELP,REFE2N,REFE2P

     Logical :: INCREP

     Integer :: mpi_rank,mpi_err,mpi_world,status(MPI_STATUS_SIZE)
     Integer :: NFIREP,tag,counter,numberConstraintsMax,exitstat
     Integer :: number_orig,task,iq,stringLength=68,size_array,mpi_string_68
     Integer, Allocatable :: buf_lambdaJob(:),buf_miuJob(:),l_temp(:,:),m_temp(:,:)

     Character(Len=6) :: string_number
     Character(Len=68) :: FILREP
     Character(Len=150) :: commande

     Real, Allocatable :: buf_qValueJob(:),q_temp(:,:)

     Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err)
     Call mpi_comm_size(MPI_COMM_WORLD,mpi_world,mpi_err)

     ! Define MPI type for a character string of length stringLength=68 and commit it
     Call mpi_type_contiguous(stringLength,MPI_CHARACTER,mpi_string_68,mpi_err)
     Call mpi_type_commit(mpi_string_68,mpi_err)

     ! Process 0 reads list of files to use as restart
     If(mpi_rank.Eq.0) Then
        numberQ=1; Call read_fileList()
        If(debug.Ge.1 .And. mpi_rank.Eq.0) Write(6,'("From hfodd_files.d - numberQ = ",i10)') numberQ
     End If

     ! Broadcasting size of collective space, and allocate space for list of filenames (at this
     ! point, only process 0 has done the allocation)
     Call MPI_Bcast(numberQ,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err)
     If(mpi_rank.Ne.0) Then
        Allocate(fichier(1:numberQ),index_file(1:numberQ))
     End If

     Call mpi_barrier(MPI_COMM_WORLD,mpi_err)

     ! Broadcasting list of files
     Call MPI_Bcast(fichier,numberQ,mpi_string_68,0,MPI_COMM_WORLD,mpi_err)
     Call MPI_Bcast(index_file,numberQ,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err)

     ! Allocate space for THO files if needed. Assume there is a directory ./tho/
     ! with all required files (this is the consequence of being in restart mode)
     If(IF_THO.Ge.1) Then
     	Allocate(fichier_tho(1:numberQ))
     	Do counter=1,numberQ
     	   Write(string_number,'(i6.6)') index_file(counter)
     	   fichier_tho(counter)='./tho/'//string_number//'.hel'
     	End Do
     End If

     ! Make a backup copy of the requested file somewhere else
     Do counter=1,numberQ
        If(Mod(mpi_rank,mpi_world).Eq.Mod(counter-1,mpi_world)) Then
           commande = 'cp '//Trim(fichier(counter))//' .'; exitstat=0
           Call system(commande)
           If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," counter = ",i10," - Executed command: ",a)') &
                                                  mpi_rank,counter,commande
        End If
     End Do

     Call mpi_barrier(MPI_COMM_WORLD,mpi_err)

     ! Reading files and obtaining characteristics of restart jobs. MPI tags are used to synchronize
     ! the order of send and receives
     If(mpi_rank.Eq.0) Then

        Allocate(qValueJob(1:numberQ,1:numberConstraints))
        Allocate(lambdaJob(1:numberQ,1:numberConstraints),miuJob(1:numberQ,1:numberConstraints))
        qValueJob(:,:)=0.0D0; lambdaJob(:,:)=0; miuJob(:,:)=0

        ! Receive the multipole constraints read from each file
        Do task=1,numberQ
           tag=task
           ! Allocate space for current buffers
           Allocate(buf_lambdaJob(numberConstraints),&
                    buf_miuJob(numberConstraints),   &
                    buf_qValueJob(numberConstraints))
           If(debug.Ge.3) Write(6,'("mpi_rank ",i10," has allocated buf_[]")') mpi_rank
           ! Receive local buffers
           tag=  numberQ+task; Call mpi_recv(buf_lambdaJob,numberConstraints,MPI_INTEGER,MPI_ANY_SOURCE,tag, &
                                                                             MPI_COMM_WORLD,status,mpi_err)
           tag=2*numberQ+task; Call mpi_recv(buf_miuJob,numberConstraints,MPI_INTEGER,MPI_ANY_SOURCE,tag, &
                                                                          MPI_COMM_WORLD,status,mpi_err)
           tag=3*numberQ+task; Call mpi_recv(buf_qValueJob,numberConstraints,MPI_DOUBLE_PRECISION,MPI_ANY_SOURCE, &
                                                                         tag,MPI_COMM_WORLD,status,mpi_err)
           If(debug.Ge.3) Write(6,'("mpi_rank ",i10," has received buf_[]")') mpi_rank
           lambdaJob(task,1:numberConstraints)=buf_lambdaJob(1:numberConstraints)
           miuJob(task,1:numberConstraints)=buf_miuJob(1:numberConstraints)
           qValueJob(task,1:numberConstraints)=buf_qValueJob(1:numberConstraints)
           Deallocate(buf_lambdaJob,buf_miuJob,buf_qValueJob)
        End Do

     Else

        Do task=1,numberQ
           If(Mod(mpi_rank,mpi_world).Eq.Mod(task-1,mpi_world)) Then
              tag=task
              FILREP=fichier(task); NFIREP=42; INCREP=.False.
              If(debug.Ge.3) Write(6,'("mpi_rank ",i10,", task ",i10,", index_file = ",i10)') &
                                   mpi_rank,task,index_file(task)
              ! Read each file
              Call read_hfodd(NFIREP,FILREP,NUMITE,NXHERM,NYHERM,NZHERM,IPCONT, &
                              ILCONT,IACONT,IMCONT,IRCONT,IRENMA,IRENIN,LIPKIN, &
                              LIPKIP,REFERN,REFERP,REDELN,REDELP,REFE2N,REFE2P, &
                              IFIBLN,INIBLN,IFIBLP,INIBLP,NMUMAX,NSIMAX,NMUCON, &
                                                                 ISHIFT,IERROR)
              ! Make sure to close file if not done before
              Inquire(unit=NFIREP,opened=INCREP); If(INCREP) Close(NFIREP)
              ! Allocate space for current buffers
              Allocate(buf_lambdaJob(numberConstraints),&
                       buf_miuJob(numberConstraints),   &
                       buf_qValueJob(numberConstraints))
              ! For each constraint iq, flex_lambda(iq) contains the lambda value, flex_miu(iq) the miu
              ! value and flex_qlm(iq) the q(lambda,miu) value
              buf_lambdaJob(:)=flex_lambda(:)
              buf_miuJob(:)=flex_miu(:)
              buf_qValueJob(:)=flex_qlm(:)
              tag=  numberQ+task; Call mpi_send(buf_lambdaJob,numberConstraints,MPI_INTEGER,0,tag, &
                                                                                MPI_COMM_WORLD,mpi_err)
              tag=2*numberQ+task; Call mpi_send(buf_miuJob,numberConstraints,MPI_INTEGER,0,tag, &
                                                                             MPI_COMM_WORLD,mpi_err)
              tag=3*numberQ+task; Call mpi_send(buf_qValueJob,numberConstraints,MPI_DOUBLE_PRECISION,0, &
                                                                            tag,MPI_COMM_WORLD,mpi_err)
              Deallocate(buf_lambdaJob,buf_miuJob,buf_qValueJob)
              If(debug.Ge.3) Write(6,'("mpi_rank ",i10," has sent buf_[]")') mpi_rank
           End If
        End Do

     End If
     If(debug.Ge.3) Write(6,'("mpi_rank = ",i10," has gone over the main loop")') mpi_rank

     If(.Not.Allocated(qValueJob)) Then
        Allocate(qValueJob(1:numberQ,1:numberConstraints))
        Allocate(lambdaJob(1:numberQ,1:numberConstraints),miuJob(1:numberQ,1:numberConstraints))
     End If
     If(debug.Ge.1) Write(6,'("mpi_rank = ",i10," has numberConstraints = ",i10)') &
                               mpi_rank, numberConstraints

     Call mpi_barrier(MPI_COMM_WORLD,mpi_err)

     ! Broadcast the full set of multipole constraints to all. First flatten 2D arrays into large vectors
     size_array=numberQ*numberConstraints
     Allocate(buf_lambdaJob(1:size_array),buf_miuJob(1:size_array),buf_qValueJob(1:size_array))
     If(mpi_rank.Eq.0) Then
        counter = 0
        Do task=1,numberQ
           Do iq=1,numberConstraints
              counter=counter+1
              buf_lambdaJob(counter)=lambdaJob(task,iq)
              buf_miuJob(counter)=miuJob(task,iq)
              buf_qValueJob(counter)=qValueJob(task,iq)
           End Do
        End Do
     Else
        buf_lambdaJob(:)=0; buf_miuJob(:)=0; buf_qValueJob(:)=0.0D0
     End If
     Call MPI_Bcast(buf_lambdaJob,size_array,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err)
     Call MPI_Bcast(buf_miuJob,size_array,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err)
     Call MPI_Bcast(buf_qValueJob,size_array,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpi_err)
     ! Then reconstruct 2D arrays from broadcast vectors
     counter = 0
     Do task=1,numberQ
        Do iq=1,numberConstraints
           counter=counter+1
           lambdaJob(task,iq)=buf_lambdaJob(counter)
           miuJob(task,iq)=buf_miuJob(counter)
           qValueJob(task,iq)=buf_qValueJob(counter)
        End Do
     End Do
     Deallocate(buf_lambdaJob,buf_miuJob,buf_qValueJob)
     If(debug.Ge.4) Then
        If(mpi_rank.Eq.0) Then
           Do task=1,numberQ
              Write(6,'("File number",i10," is ",a68)') task,fichier(task)
              Do iq=1,numberConstraints
                 Write(6,'("l = ",i4," m = ",i4," qlm = ",f10.5)') &
                       lambdaJob(task,iq),miuJob(task,iq),qValueJob(task,iq)
              End Do
           End Do
         End If
     End If
     ! Set original number of constraints and original number of "grid" points
     numberQ_orig = numberQ
     numberConstraints_orig = numberConstraints
     neck_numberPoints = 1

     ! Releasing derived MPI datatype
     Call mpi_type_free(mpi_string_68,mpi_err)

   End Subroutine define_restartPointsQlm

   !---------------------------------------------------------------------!
   !  Function returning the job number embedded in a HFODD record file  !
   !  name, used when modeMPI_deformation = 3 only                       !
   !---------------------------------------------------------------------!
   Integer Function extract_index(filename)

     Character(Len=68), INTENT(IN) :: filename
     Integer :: i

     If(filename(1:1).Eq.'H') Then
        Read(filename,'("HFODD_",i8.8,".REC")') i
        extract_index = i
     End If

   End Function extract_index

   !---------------------------------------------------------------------!
   !  Subroutine reading the list of files to restart from in the case   !
   !  where modeMPI_deformation = 3. Executed only from proc 0           !
   !  Outputs: - numberQ, total number of points on the mesh             !
   !           - fichier(:), list of files to restart from               !
   !---------------------------------------------------------------------!
   Subroutine read_fileList()

     Integer :: unit,counter
     Character(Len=68) :: filename

     unit=41

     ! First pass through the file to count the number of files: defines the initial "mesh"
     counter=0; filename=''
     Open(unit,FILE='hfodd_files.d',STATUS='UNKNOWN',FORM='FORMATTED')
     Do While(trim(filename) .Ne. 'ALL_DONE' .And. trim(filename) .Ne. 'ALL-DONE')
        Read(unit,*) filename
        counter=counter+1
     End Do
     Close(unit)

     numberQ=counter-1

     ! Second pass to actually read the filename
     filename=''
     Allocate(fichier(1:numberQ),index_file(1:numberQ))
     Open(unit,FILE='hfodd_files.d',STATUS='UNKNOWN',FORM='FORMATTED')
     Do counter=1,numberQ
        Read(unit,*) filename
        index_file(counter)=extract_index(filename)
        fichier(counter)='./restart/'//filename(1:60)
     End Do
     Close(unit)

   End Subroutine read_fileList

   !---------------------------------------------------------------------!
   !  Subroutine mapping the default deformation grid to the restart     !
   !  deformation grid. Only used when modeMPI_deformation = 2. For this !
   !  mapping to function, it is assumed that the restart collective     !
   !  space is equal or larger than the original one: there may be more  !
   !  constraints, but not less; the number of points for the shared     !
   !  collective variable may be different.                              !
   !  Outputs: - ii=mapping_qlm(jj) gives the index ii of the original   !
   !             grid used to restart the calculation on point jj of the !
   !             restart grid                                            !
   !           - numberQ, numberConstraints are overwritten and contain  !
   !             the characteristics of the restart collective space     !
   !           - qValueJob(:), lambdaJob(:) and miuJob(:) are similarly  !
   !             associated with the restart grid                        !
   !  TODO: THE CURRENT METHOD MAY BE TAXING THE OS WHEN THE NEW GRID IS !
   !        LARGER THAN THE OLD ONE. INDEED, THE SAME FILE MAY BE READ   !
   !	    BY MANY DIFFERENT PROCESSES, WHICH IS NOT EFFICIENT.         !
   !---------------------------------------------------------------------!
   Subroutine mapping_gridQlm(IF_THO,LIPKIN,LIPKIP)

     Integer, INTENT(IN) :: IF_THO,LIPKIN,LIPKIP
     Integer :: ii,jj,index_optimal,iq,jq,shift,lambda_restart,miu_restart, &
                lambda,miu,mpi_rank,mpi_err,mpi_world,exitstat,pass_number
     Real :: distance_min,distance,qlm_restart,qlm
     Character(Len=132) :: directory_name
     Character(Len=68) :: fichier_courant,fichier_optimal
     Character(Len=150) :: commande

     Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err)
     Call mpi_comm_size(MPI_COMM_WORLD,mpi_world,mpi_err)

     ! Default mapping: all new calculations restart from HFODD_00000001.REC
     Allocate(mapping_qlm(1:numberQ_restart))
     Do jj=1,numberQ_restart
        mapping_qlm(jj)=1
     End Do

     ! Move all restart files into current directory
     If(mpi_rank.Eq.0.And.modeMPI_neck.Ne.2) Then
        commande = 'mv restart restart_old'; exitstat=0
        Call system(commande)
        If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') &
                                               mpi_rank,commande
        commande = 'mkdir restart'; exitstat=0
        Call system(commande)
        If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') &
                                               mpi_rank,commande
        If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then
           commande = 'mv lic lic_old'; exitstat=0
           Call system(commande)
           If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') &
                                                  mpi_rank,commande
           commande = 'mkdir lic'; exitstat=0
           Call system(commande)
           If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') &
                                                  mpi_rank,commande
        End If
        ! In batch mode, restart files will later be overwritten with record files. We make an
        ! additional copy here
        If(batch_mode.Eq.1) Then
           commande = 'mkdir restart_tmp'; exitstat=0
           Call system(commande)
           If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') &
                                                  mpi_rank,commande
           If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then
              commande = 'mkdir lic_tmp'; exitstat=0
              Call system(commande)
              If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') &
                                                     mpi_rank,commande
           End If
        End If
     End If

     ! Wait here for all processes
     Call mpi_barrier(MPI_COMM_WORLD,mpi_err)

     Allocate(fichier_rec_old(1:numberQ_restart),fichier_rec_new(1:numberQ_restart))
     If(IF_THO.Ge.1) Allocate(fichier_tho_old(1:numberQ_restart),fichier_tho_new(1:numberQ_restart))
     If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Allocate(fichier_lic_old(1:numberQ_restart),fichier_lic_new(1:numberQ_restart))
     If(batch_mode.Eq.1) Then
        Allocate(fichier_rec_tmp(1:numberQ_restart))
        If(IF_THO.Ge.1) Allocate(fichier_tho_tmp(1:numberQ_restart))
        If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Allocate(fichier_lic_tmp(1:numberQ_restart))
     End If

     ! Loop over the new grid
     Do jj=1,numberQ_restart

        index_optimal=1; distance_min=1.D9

        ! Loop over the old grid
        Do ii=1,numberQ

           distance = 0.0D0

           ! Loop over the old constraints
           Do iq=1,numberConstraints

              ! Both old and new constraints were read using linked lists. If Nold is the
              ! number of old constraints and Nnew that of the new constraints, we have
              ! Nnew >= Nold (necessarily) and the shared Nold constraints are stored at
              ! the end of lambdaJob(:), lambdaJob_restart(:), etc. To ge the proper
              ! correspondence, we use the shift below
              shift = iq + (numberConstraints_restart - numberConstraints)
              lambda_restart=lambdaJob_restart(jj,shift)
              miu_restart= miuJob_restart(jj,shift)
              qlm_restart=qValueJob_restart(jj,shift)

              lambda=lambdaJob(ii,iq)
              miu=miuJob(ii,iq)
              qlm=qValueJob(ii,iq)

              ! Measuring the distance between point ii (old) and jj (new)
              ! in terms of multipole moments
              If(lambda.Eq.lambda_restart .And. miu.Eq.miu_restart) Then
                 distance = distance + (qlm - qlm_restart)**2
              End If

           End Do

           ! Looking for the point index_optimal that minimizes the distance
           If(distance.Le.distance_min) Then
              distance_min=distance; index_optimal = ii
           End If

        End Do

        ! mv proper file into restart directory with proper name
        If(modeMPI_neck.Ne.2)Then
           Write(fichier_optimal,'("./restart_old/HFODD_",I8.8,".REC")') index_optimal
           Write(fichier_courant,'("./restart/HFODD_",I8.8,".REC")') jj
           fichier_rec_old(jj)=Trim(fichier_optimal)
           fichier_rec_new(jj)=Trim(fichier_courant)
           If(IF_THO.Ge.1) Then
              Write(fichier_optimal,'("./tho/t",I6.6,".hel")') index_optimal
              Write(fichier_courant,'("./t",I6.6,".hel")') jj
              fichier_tho_old(jj)=Trim(fichier_optimal)
              fichier_tho_new(jj)=Trim(fichier_courant)
           End If
           If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then
              Write(fichier_optimal,'("./lic_old/HFODD_",I8.8,".LIC")') index_optimal
              Write(fichier_courant,'("./lic/HFODD_",I8.8,".LIC")') jj
              fichier_lic_old(jj)=Trim(fichier_optimal)
              fichier_lic_new(jj)=Trim(fichier_courant)
           End If
           mapping_qlm(jj)=jj
           ! In batch mode, we need to make an extra copy of all these restart files, since they will
           ! be overwritten later
           If(batch_mode.Eq.1) Then
              Write(fichier_courant,'("./restart_tmp/HFODD_",I8.8,".REC")') jj
              fichier_rec_tmp(jj)=Trim(fichier_courant)
              If(IF_THO.Eq.1) Then
                 Write(fichier_courant,'("./t",I6.6,".hel")') jj
                 fichier_tho_tmp(jj)=Trim(fichier_courant)
              End If
              If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then
                 Write(fichier_courant,'("./lic_tmp/HFODD_",I8.8,".LIC")') jj
                 fichier_lic_tmp(jj)=Trim(fichier_courant)
              End If
           End If
        Else
           mapping_qlm(jj)=index_optimal
        End If

        If(mpi_rank.Eq.0.And.debug.Ge.4) Then
           Write(6,'("mpi_rank = ",i10," mapping_qlm - jj = ",i10," index_optimal = ",i10)') &
                     mpi_rank,jj,index_optimal
           Do jq=1,numberConstraints_restart
              Write(6,'("RESTART  - mpi_rank = ",i10," l=",i2,"m=",i2," Qlm=",f10.5)') &
                                    mpi_rank,lambdaJob_restart(jj,jq),                 &
                                                miuJob_restart(jj,jq),                 &
                                             qValueJob_restart(jj,jq)
           End do
           Do iq=1,numberConstraints
              Write(6,'("ORIGINAL - mpi_rank = ",i10," l=",i2,"m=",i2," Qlm=",f10.5)') &
                                    mpi_rank,lambdaJob(index_optimal,iq),              &
                                                miuJob(index_optimal,iq),              &
                                 qValueJob(index_optimal,iq)
           End do
        End If

     End Do

     If(modeMPI_neck.Ne.2) Then

     	! Loop over the new grid
     	Do jj=1,numberQ_restart
     	   If(Mod(mpi_rank,mpi_world).Eq.Mod(jj-1,mpi_world)) Then
     	      commande = 'cp '//fichier_rec_old(jj)//' '//fichier_rec_new(jj); exitstat=0
     	      Call system(commande)
     	      If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
     	                                             mpi_rank,jj,commande
     	      If(IF_THO.Ge.1) Then
     	         commande = 'cp '//fichier_tho_old(jj)//' '//fichier_tho_new(jj); exitstat=0
     	         Call system(commande)
     	         If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
     	                                                mpi_rank,jj,commande
     	      End If
     	      If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then
     	         commande = 'cp '//fichier_lic_old(jj)//' '//fichier_lic_new(jj); exitstat=0
     	         Call system(commande)
     	         If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
     	                                                mpi_rank,jj,commande
     	      End If
     	      ! Make extra copy in batch mode
     	      If(batch_mode.Eq.1) Then
                 commande = 'cp '//fichier_rec_old(jj)//' '//fichier_rec_tmp(jj); exitstat=0
                 Call system(commande)
                 If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
                                                       mpi_rank,jj,commande
                 If(IF_THO.Ge.1) Then
                    commande = 'cp '//fichier_tho_old(jj)//' '//fichier_tho_tmp(jj); exitstat=0
                    Call system(commande)
                    If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
                                                          mpi_rank,jj,commande
                 End If
     	         If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then
     	            commande = 'cp '//fichier_lic_old(jj)//' '//fichier_lic_tmp(jj); exitstat=0
     	            Call system(commande)
     	            If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
     	                                                   mpi_rank,jj,commande
     	         End If
     	      End If
     	   End If
     	End Do

        ! Wait here for all processes
        Call mpi_barrier(MPI_COMM_WORLD,mpi_err)

        ! Check that all files have truly been copied
        Allocate(table_fichiers(1:numberQ_restart)); table_fichiers(:)=0
        directory_name='restart'

        pass_number=1; filled=.False.
        Do While (.Not.filled)

           If(mpi_rank.Eq.0.And.debug.Ge.2) Write(6,'("Pass=",i5)') pass_number
           Call check_directory(directory_name,numberQ_restart)

           Call fill_holes(numberQ_restart,IF_THO,LIPKIN,LIPKIP)

           pass_number=pass_number+1
           If(pass_number.Gt.4) Exit ! avoid infinite loops

        End Do
        Deallocate(table_fichiers)

        ! Check that all files in restart_tmp/ have truly been copied (batch mode)
        If(batch_mode.Eq.1) Then

           Allocate(table_fichiers(1:numberQ_restart)); table_fichiers(:)=0
           directory_name='restart_tmp'

           pass_number=1; filled=.False.
           Do While (.Not.filled)

              If(mpi_rank.Eq.0.And.debug.Ge.2) Write(6,'("Pass =",i5)') pass_number
              Call check_directory(directory_name,numberQ_restart)

              Call fill_holes(numberQ_restart,IF_THO,LIPKIN,LIPKIP)

              pass_number=pass_number+1
              If(pass_number.Gt.4) Exit ! avoid infinite loops

           End Do
           Deallocate(table_fichiers)

     End If

     End If ! End modeMPI_neck.Ne.2

     Deallocate(fichier_rec_old,fichier_rec_new)
     If(IF_THO.Ge.1) Deallocate(fichier_tho_old,fichier_tho_new)
     If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Deallocate(fichier_lic_old,fichier_lic_new)
     If(batch_mode.Eq.1) Then
        Deallocate(fichier_rec_tmp)
        If(IF_THO.Ge.1) Deallocate(fichier_tho_tmp)
        If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Deallocate(fichier_lic_tmp)
     End If

     Deallocate(qValueJob,lambdaJob,miuJob)

     numberQ=numberQ_restart
     numberConstraints=numberConstraints_restart

     Allocate(qValueJob(1:numberQ,1:numberConstraints))
     Allocate(lambdaJob(1:numberQ,1:numberConstraints),miuJob(1:numberQ,1:numberConstraints))
     qValueJob(:,:)=qValueJob_restart(:,:)
     lambdaJob(:,:)=lambdaJob_restart(:,:); miuJob(:,:)=miuJob_restart(:,:)
     Deallocate(qValueJob_restart,lambdaJob_restart,miuJob_restart)

   End Subroutine mapping_gridQlm

   !---------------------------------------------------------------------!
   !  Subroutine mapping the default deformation grid to the restart     !
   !  deformation grid. Only used when modeMPI_deformation >= 3. In this !
   !  case, the original collective space is arbitrary compared to the   !
   !  restart one: it may have more or fewer collective variables, the   !
   !  mesh can be irregular, etc. This is an option meant to provide     !
   !  flexibility.                                                       !
   !  Outputs: - ii=mapping_qlm_new(jj) gives as before the index ii of  !
   !             the original grid used to restart the calculation on    !
   !             point jj of the restart grid                            !
   !           - numberQ, numberConstraints are overwritten and contain  !
   !             the characteristics of the restart collective space     !
   !           - qValueJob(:), lambdaJob(:) and miuJob(:) are similarly  !
   !             associated with the restart grid                        !
   !---------------------------------------------------------------------!
   Subroutine mapping_gridQlm_new(IF_THO,LIPKIN,LIPKIP)

     Integer, INTENT(IN) :: IF_THO,LIPKIN,LIPKIP
     Integer :: ii,jj,index_optimal,iq,jq,lambda_restart,miu_restart, &
                lambda,miu,active,mpi_rank,mpi_err,mpi_world,exitstat,&
                pass_number
     Real :: distance_min,distance,qlm_restart,qlm
     Character(Len=132) :: directory_name
     Character(Len=68) :: fichier_courant,fichier_optimal
     Character(Len=150) :: commande

     Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err)
     Call mpi_comm_size(MPI_COMM_WORLD,mpi_world,mpi_err)

     ! Default mapping: all new calculations restart from HFODD_00000001.REC
     Allocate(mapping_qlm_new(1:numberQ_restart))
     Do jj=1,numberQ_restart
        mapping_qlm_new(jj)=1
     End Do

     ! Move all restart files into current directory
     If(mpi_rank.Eq.0) Then
        commande = 'mv restart restart_old'; exitstat=0
        Call system(commande)
        If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') &
                                               mpi_rank,commande
        commande = 'mkdir restart'; exitstat=0
        Call system(commande)
        If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') &
                                               mpi_rank,commande
        If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then
           commande = 'mv lic lic_old'; exitstat=0
           Call system(commande)
           If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') &
                                                  mpi_rank,commande
           commande = 'mkdir lic'; exitstat=0
           Call system(commande)
           If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') &
                                                  mpi_rank,commande
        End If
        ! In batch mode, restart files will later be overwritten with record files. We make an
        ! additional copy here
        If(batch_mode.Eq.1) Then
           commande = 'mkdir restart_tmp'; exitstat=0
           Call system(commande)
           If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') &
                                                  mpi_rank,commande
           If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then
              commande = 'mkdir lic_tmp'; exitstat=0
              Call system(commande)
              If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') &
                                                     mpi_rank,commande
           End If
        End If
     End If

     ! Wait here for all processes
     Call mpi_barrier(MPI_COMM_WORLD,mpi_err)

     ! Allocate filenames
     Allocate(fichier_rec_old(1:numberQ_restart),fichier_rec_new(1:numberQ_restart))
     If(IF_THO.Ge.1) Allocate(fichier_tho_old(1:numberQ_restart),fichier_tho_new(1:numberQ_restart))
     If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Allocate(fichier_lic_old(1:numberQ_restart),fichier_lic_new(1:numberQ_restart))
     If(batch_mode.Eq.1) Then
        Allocate(fichier_rec_tmp(1:numberQ_restart))
        If(IF_THO.Ge.1) Allocate(fichier_tho_tmp(1:numberQ_restart))
        If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Allocate(fichier_lic_tmp(1:numberQ_restart))
     End If

     ! Loop over the new grid
     Do jj=1,numberQ_restart

        index_optimal=1; distance_min=1.D9

        ! Loop over the old grid
        Do ii=1,numberQ

           distance=0.0D0; active=0

           ! Loop over the new collective variables (=constraints)
           Do jq=1,numberConstraints_restart
              lambda_restart=lambdaJob_restart(jj,jq)
              miu_restart=miuJob_restart(jj,jq)
              qlm_restart=qValueJob_restart(jj,jq)
              ! Loop over the old collective variables (=constraints)
              Do iq=1,numberConstraints
                 lambda=lambdaJob(ii,iq)
                 miu=miuJob(ii,iq)
                 qlm=qValueJob(ii,iq)
                 ! Measuring the distance between point ii (old) and jj (new)
                 ! in terms of multipole moments
                 If(lambda.Eq.lambda_restart .And. miu.Eq.miu_restart) Then
                    active=active+1
                    distance = distance + (qlm - qlm_restart)**2
                 End If
              End Do
           End Do

           ! Looking for the point index_optimal that minimizes the distance
           If(distance.Le.distance_min.And.active.Gt.0) Then
              distance_min=distance; index_optimal = ii
           End If

        End Do

        ! move proper file into restart directory with proper name
        Write(fichier_optimal,'("./restart_old/HFODD_",I8.8,".REC")') index_file(index_optimal)
        Write(fichier_courant,'("./restart/HFODD_",I8.8,".REC")') jj
        fichier_rec_old(jj)=Trim(fichier_optimal)
        fichier_rec_new(jj)=Trim(fichier_courant)
        If(IF_THO.Ge.1) Then
           Write(fichier_optimal,'("./tho/t",I6.6,".hel")') index_file(index_optimal)
           Write(fichier_courant,'("./t",I6.6,".hel")') jj
           fichier_tho_old(jj)=Trim(fichier_optimal)
           fichier_tho_new(jj)=Trim(fichier_courant)
        End If
        If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then
           Write(fichier_optimal,'("./lic_old/HFODD_",I8.8,".LIC")') index_optimal
           Write(fichier_courant,'("./lic/HFODD_",I8.8,".LIC")') jj
           fichier_lic_old(jj)=Trim(fichier_optimal)
           fichier_lic_new(jj)=Trim(fichier_courant)
        End If
        ! In batch mode, we need to make an extra copy of all these restart files, since they will
        ! be overwritten later
        If(batch_mode.Eq.1) Then
           Write(fichier_courant,'("./restart_tmp/HFODD_",I8.8,".REC")') jj
           fichier_rec_tmp(jj)=Trim(fichier_courant)
           If(IF_THO.Eq.1) Then
              Write(fichier_courant,'("./t",I6.6,".hel")') jj
              fichier_tho_tmp(jj)=Trim(fichier_courant)
           End If
           If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then
              Write(fichier_courant,'("./lic_tmp/HFODD_",I8.8,".LIC")') jj
              fichier_lic_tmp(jj)=Trim(fichier_courant)
           End If
        End If
        mapping_qlm_new(jj)=jj

        If(mpi_rank.Eq.0.And.debug.Ge.4) Then
           Write(6,'("mpi_rank = ",i10," mapping_qlm_new - jj = ",i10," index_optimal = ",i10)') &
                     mpi_rank,jj,index_file(index_optimal)
           Do jq=1,numberConstraints_restart
              Write(6,'("RESTART  - mpi_rank = ",i10," l=",i2,"m=",i2," Qlm=",f10.5)') &
                                    mpi_rank,lambdaJob_restart(jj,jq),                 &
                                             miuJob_restart(jj,jq),                    &
                                 qValueJob_restart(jj,jq)
           End do
           Do iq=1,numberConstraints
              Write(6,'("ORIGINAL - mpi_rank = ",i10," l=",i2,"m=",i2," Qlm=",f10.5)') &
                                    mpi_rank,lambdaJob(index_optimal,iq),              &
                                              miuJob(index_optimal,iq),                &
                                 qValueJob(index_optimal,iq)
           End do
        End If

     End Do

     ! Loop over the new grid
     Do jj=1,numberQ_restart
        If(Mod(mpi_rank,mpi_world).Eq.Mod(jj-1,mpi_world)) Then
           commande = 'cp '//fichier_rec_old(jj)//' '//fichier_rec_new(jj); exitstat=0
           Call system(commande)
           If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
                                                 mpi_rank,jj,commande
           If(IF_THO.Ge.1) Then
              commande = 'cp '//fichier_tho_old(jj)//' '//fichier_tho_new(jj); exitstat=0
              Call system(commande)
              If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
                                                    mpi_rank,jj,commande
           End If
     	   If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then
     	      commande = 'cp '//fichier_lic_old(jj)//' '//fichier_lic_new(jj); exitstat=0
     	      Call system(commande)
     	      If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
     	                                             mpi_rank,jj,commande
     	   End If
     	   ! Make extra copy in batch mode
     	   If(batch_mode.Eq.1) Then
              commande = 'cp '//fichier_rec_old(jj)//' '//fichier_rec_tmp(jj); exitstat=0
              Call system(commande)
              If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
                                                    mpi_rank,jj,commande
              If(IF_THO.Ge.1) Then
                 commande = 'cp '//fichier_tho_old(jj)//' '//fichier_tho_tmp(jj); exitstat=0
                 Call system(commande)
                 If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
                                                       mpi_rank,jj,commande
              End If
     	      If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then
     	         commande = 'cp '//fichier_lic_old(jj)//' '//fichier_lic_tmp(jj); exitstat=0
     	         Call system(commande)
     	         If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
     	                                                mpi_rank,jj,commande
     	      End If
     	   End If
        End If
     End Do

     ! Wait here for all processes
     Call mpi_barrier(MPI_COMM_WORLD,mpi_err)

     ! Check that all files in restart/ have truly been copied
     Allocate(table_fichiers(1:numberQ_restart)); table_fichiers(:)=0
     directory_name='restart'

     pass_number=1; filled=.False.
     Do While (.Not.filled)

        If(mpi_rank.Eq.0.And.debug.Ge.2) Write(6,'("Pass =",i5)') pass_number
        Call check_directory(directory_name,numberQ_restart)

        Call fill_holes(numberQ_restart,IF_THO,LIPKIN,LIPKIP)

        pass_number=pass_number+1
        If(pass_number.Gt.4) Exit ! avoid infinite loops

     End Do
     Deallocate(table_fichiers)

     ! Check that all files in restart_tmp/ have truly been copied (batch mode)
     If(batch_mode.Eq.1) Then

        Allocate(table_fichiers(1:numberQ_restart)); table_fichiers(:)=0
        directory_name='restart_tmp'

        pass_number=1; filled=.False.
        Do While (.Not.filled)

           If(mpi_rank.Eq.0.And.debug.Ge.2) Write(6,'("Pass =",i5)') pass_number
           Call check_directory(directory_name,numberQ_restart)

           Call fill_holes(numberQ_restart,IF_THO,LIPKIN,LIPKIP)

           pass_number=pass_number+1
           If(pass_number.Gt.4) Exit ! avoid infinite loops

        End Do
        Deallocate(table_fichiers)
     End If

     Deallocate(fichier_rec_old,fichier_rec_new)
     If(IF_THO.Ge.1) Deallocate(fichier_tho_old,fichier_tho_new)
     If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Deallocate(fichier_lic_old,fichier_lic_new)
     If(batch_mode.Eq.1) Then
        Deallocate(fichier_rec_tmp)
        If(IF_THO.Ge.1) Deallocate(fichier_tho_tmp)
        If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Deallocate(fichier_lic_tmp)
     End If

     Deallocate(qValueJob,lambdaJob,miuJob)

     numberQ=numberQ_restart
     numberConstraints=numberConstraints_restart
     If(mpi_rank.Eq.0.And.debug.Ge.2) &
        Write(6,'("In mapping_gridQlm_new - numberQ =",i10," numberConstraints = ",i5)') numberQ, numberConstraints

     Allocate(qValueJob(1:numberQ,1:numberConstraints))
     Allocate(lambdaJob(1:numberQ,1:numberConstraints),miuJob(1:numberQ,1:numberConstraints))
     qValueJob(:,:)=qValueJob_restart(:,:)
     lambdaJob(:,:)=lambdaJob_restart(:,:); miuJob(:,:)=miuJob_restart(:,:)
     Deallocate(qValueJob_restart,lambdaJob_restart,miuJob_restart)

   End Subroutine mapping_gridQlm_new

   !---------------------------------------------------------------------!
   !  Subroutine mapping the default neck grid to the restart neck grid. !
   !  Outputs: - ii=mapping_neck(jj) gives the index ii of the file used !
   !             to restart the calculation of point jj on the restart   !
   !             grid.                                                   !
   !           - neck_value(:), neck_numberPoints are overwritten and    !
   !             now correspond to the restart grid                      !
   !---------------------------------------------------------------------!
   Subroutine mapping_gridNeck()

     Integer :: ii,jj,index_optimal
     Real :: qn,qn_restart,distance,distance_min

     ! Default mapping: all new calculations restart from HFODD_00000001.REC
     Allocate(mapping_neck(1:neck_numberPoints_restart))
     Do jj=1,neck_numberPoints_restart
        mapping_neck(jj)=1
     End Do

     Do jj=1,neck_numberPoints_restart

        index_optimal=1; distance_min=1.0D9

        Do ii=1,neck_numberPoints

           qn_restart = neck_value_restart(jj)
           qn = neck_value(ii)

           distance = (qn - qn_restart)**2

           If(distance.Le.distance_min) Then
              distance_min=distance; index_optimal = ii
           End If

        End Do

        mapping_neck(jj)=index_optimal
        If(debug.Ge.1) Write(6,'("mapping_gridNeck - jj = ",i10," index_optimal = ",i10)') &
                                                     jj,index_optimal

     End Do

     neck_numberPoints = neck_numberPoints_restart
     Deallocate(neck_value); Allocate(neck_value(1:neck_numberPoints))
     neck_value(:) = neck_value_restart(:)
     Deallocate(neck_value_restart)

   End Subroutine mapping_gridNeck

   !---------------------------------------------------------------------!
   !  Subroutine defining the tasks                                      !
   !  Possible combinations:                                             !
   !     modeMPI_deformation  1  1  2  2  2  3                           !
   !     modeMPI_neck         0  1  0  1  2  0                           !
   !---------------------------------------------------------------------!
   Subroutine mpi_constructJobList(NUMITE,NXHERM,NYHERM,NZHERM,        &
                                   IPCONT,ILCONT,IACONT,IMCONT,IRCONT, &
                                                        IRENMA,IRENIN, &
                                                        LIPKIN,LIPKIP, &
                                                 REFERN,REFERP,REDELN, &
                                                 REDELP,REFE2N,REFE2P, &
                                          IFIBLN,INIBLN,IFIBLP,INIBLP, &
                                          NMUMAX,NSIMAX,NMUCON,ISHIFT, &
                                                               IF_THO, &
                                                               IERROR)

     Integer, INTENT(IN) :: NUMITE,NXHERM,NYHERM,NZHERM,IPCONT,ILCONT,IACONT, &
                            IMCONT,IRCONT,IRENMA,IRENIN,LIPKIN,LIPKIP,IFIBLN, &
                            INIBLN,IFIBLP,INIBLP,NMUMAX,NSIMAX,NMUCON,ISHIFT, &
                            IF_THO
     Integer, INTENT(INOUT) :: IERROR
     Real, INTENT(INOUT) :: REFERN,REFERP,REDELN,REDELP,REFE2N,REFE2P

     Integer :: mpi_rank, mpi_err
     Integer :: iz,in,iq,nCommon,nGridOrig,indexJob_orig,iCommon,force,temperature,idxNeck, &
                idxDefo,iGrid,indexJob,zNumber,nNumber,i_neck,ii,i_shel,i_stat,i_freq,i_defo
     Integer, Allocatable :: zVector(:),nVector(:),mapping_oldGrid(:,:),indexJobVector_orig(:,:)

     Real :: xzNumber, xnNumber
     Real, Allocatable :: xzVector(:),xnVector(:)

     Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err)

     If(debug.Ge.1 .And. mpi_rank .Eq.0) &
        Write(6,'("modeMPI_deformation=",i3," modeMPI_neck=",i3," modeMPI_basis=",i3)') &
                          modeMPI_deformation,modeMPI_neck,modeMPI_basis

     ! Define deformation and neck grids
     If(modeMPI_deformation.Le.2) Then
        Call define_gridQlm()
        If(modeMPI_deformation.Eq.2) Then
           Call define_gridQlm_new()
           Call mapping_gridQlm(IF_THO,LIPKIN,LIPKIP)
        End If
        ! If constraints on the size of the neck are active, we define the grids
        If(modeMPI_neck.Ge.1) Then
           Call define_gridNeck()
           If(modeMPI_neck.Eq.1) Then
              ! Default mapping if no constraint on the size of the neck is requested, neck_numberPoints
              ! is preset to 1 in routine mpi_initializeParallelData() of module hfodd_mpiio.f90
              neck_numberPoints_orig=1
              Allocate(mapping_neck(1:neck_numberPoints))
              Do ii=1,neck_numberPoints
                 mapping_neck(ii)=1
              End Do
           Else
              Call define_gridNeck_new()
              Call mapping_gridNeck()
           End If
     Else
           ! Default mapping if no constraint on the size of the neck is requested, neck_numberPoints
           ! is preset to 1 in routine mpi_initializeParallelData() of module hfodd_mpiio.f90
           neck_numberPoints_orig=1
           Allocate(mapping_neck(1:neck_numberPoints))
           Do ii=1,neck_numberPoints
              mapping_neck(ii)=1
           End Do
        End If
     Else
        Call define_restartPointsQlm(NUMITE,NXHERM,NYHERM,NZHERM,IPCONT,ILCONT,IACONT, &
                                     IMCONT,IRCONT,IRENMA,IRENIN,LIPKIN,LIPKIP,REFERN, &
                                     REFERP,REDELN,REDELP,REFE2N,REFE2P,IFIBLN,INIBLN, &
                                     IFIBLP,INIBLP,NMUMAX,NSIMAX,NMUCON,ISHIFT,IF_THO, &
                                                                               IERROR)
        If(modeMPI_deformation.Eq.3) Call define_gridQlm_new()
        If(modeMPI_deformation.Eq.4) Call define_pathQlm_new()
        Call mapping_gridQlm_new(IF_THO,LIPKIN,LIPKIP)
        ! Default mapping since no constraint on the size of the neck is allowed
        neck_numberPoints_orig=neck_numberPoints
        Allocate(mapping_neck(1:neck_numberPoints))
        Do ii=1,neck_numberPoints
           mapping_neck(ii)=ii
        End Do
     End If

     ! default allocations of mapping arrays: if no restart is required, these mappings simply
     ! provide identity vectors (v(i) = i)
     If(modeMPI_deformation.Le.1) Then
        Allocate(mapping_qlm(1:numberQ))
        Do ii=1,numberQ
           mapping_qlm(ii)=ii
           If(debug.Ge.4 .And. mpi_rank .Eq.0) Write(6,'("ii = ",i10," mapping = ",i10)') ii,mapping_qlm(ii)
        End Do
     End If

     ! Define nucleus grid: integers
     Allocate(zVector(1:proton_stepNumber),nVector(1:neutron_stepNumber))
     Do iz=1,proton_stepNumber
        zVector(iz) = proton_initValue + (iz-1)*proton_stepValue
     End Do
     Do in=1, neutron_stepNumber
        nVector(in) = neutron_initValue + (in-1)*neutron_stepValue
     End Do

     ! Define nucleus grid: float
     If(x_proton_initValue.Gt.0.0D0) Then
        Allocate(xzVector(1:proton_stepNumber))
        Do iz=1,proton_stepNumber
           xzVector(iz) = x_proton_initValue + (iz-1)*x_proton_stepValue
        End Do
     End If
     If(x_neutron_initValue.Gt.0.0D0) Then
        Allocate(xnVector(1:neutron_stepNumber))
        Do in=1, neutron_stepNumber
           xnVector(in) = x_neutron_initValue + (in-1)*x_neutron_stepValue
        End Do
     End If

     ! Default deformation mode
     If(modeMPI_deformation.Ge.1.And.modeMPI_deformation.Le.2) Then

        NOFJOB = proton_stepNumber * neutron_stepNumber * counterForce * numberQ &
               * counterTemperature * neck_numberPoints

        Allocate(IZ_JOB(1:NOFJOB),IN_JOB(1:NOFJOB),indexQjob(1:NOFJOB),indexRestartJob(1:NOFJOB),&
                 indexForce(1:NOFJOB),indexTemperature(1:NOFJOB),indexNeck(1:NOFJOB))
        Allocate(XIZJOB(1:NOFJOB),XINJOB(1:NOFJOB)); XIZJOB(:)=-1.0D0; XINJOB(:)=-1.0D0

        nCommon = proton_stepNumber * neutron_stepNumber * counterForce * counterTemperature
        nGridOrig = numberQ * neck_numberPoints_orig
        Allocate(indexJobVector_orig(1:nGridOrig, 1:nCommon))

        indexJob_orig = 0

        If(debug.Ge.1 .And. mpi_rank .Eq.0) Then
           Write(6,'("Sanity checks:")')
           Write(6,'("proton_stepNumber=",i10," neutron_stepNumber=",i10," numberConstraints=",i10)') &
                      proton_stepNumber,neutron_stepNumber,numberConstraints
           Write(6,'("counterForce=",i10," counterTemperature=",i10," restartTemperature=",i10)') &
                      counterForce,counterTemperature,restartTemperature
           Write(6,'("numberQ=",i10," numberQ_orig=",i10," neck_numberPoints=",i10," neck_numberPoints_orig=",i10)') &
                      numberQ,numberQ_orig,neck_numberPoints,neck_numberPoints_orig
        End If

        If(restartTemperature == 0) Then

           iCommon=0

           Allocate(mapping_oldGrid(1:numberQ,1:neck_numberPoints_orig))

           Do iz=1,proton_stepNumber
              Do in=1,neutron_stepNumber
                 Do force=1,counterForce
                    Do temperature=1,counterTemperature

                       iCommon=iCommon+1; iGrid=0

                       Do idxNeck=1,neck_numberPoints_orig
                          Do idxDefo=1,numberQ

                             iGrid=iGrid+1
                             mapping_oldGrid(idxDefo,idxNeck)=iGrid

                             indexJob_orig=indexJob_orig+1
                             indexJobVector_orig(iGrid,iCommon)=indexJob_orig

                          End Do
                       End Do

                    End Do
                 End Do
              End Do
           End Do

           iCommon=0; indexJob=0

           Do iz=1,proton_stepNumber
              zNumber = zVector(iz)
              If(x_proton_initValue.Gt.0.0D0) xzNumber=xzVector(iz)
              Do in=1,neutron_stepNumber
                 nNumber = nVector(in)
                 If(x_neutron_initValue.Gt.0.0D0) xnNumber=xnVector(in)
                 Do force=1,counterForce

                    iCommon = iCommon + 1

                    Do temperature=1,counterTemperature
                       Do i_neck=1,neck_numberPoints
                          Do ii=1,numberQ

                             indexJob=indexJob+1
                             indexQjob(indexJob)=ii

                             indexForce(indexJob)=force
                             indexTemperature(indexJob)=temperature
                             indexNeck(indexJob)=i_neck

                             IZ_JOB(indexJob)=zNumber
                             IN_JOB(indexJob)=nNumber

                             If(x_proton_initValue.Gt.0.0D0) Then
                                XIZJOB(indexJob)=xzNumber
                             Else
                                XIZJOB(indexJob)=Real(zNumber,Kind=Kind(1.0))
                             End If
                             If(x_neutron_initValue.Gt.0.0D0) Then
                                XINJOB(indexJob)=xnNumber
                             Else
                                XINJOB(indexJob)=Real(nNumber,Kind=Kind(1.0))
                             End If

                             idxDefo=mapping_qlm(ii)
                             idxNeck=mapping_neck(i_neck)
                             iGrid=mapping_oldGrid(idxDefo,idxNeck)

                             indexRestartJob(indexJob)=indexJobVector_orig(iGrid,iCommon)

                             If(debug.Ge.4 .And. mpi_rank .Eq.0) Then
                                Do iq=1,numberConstraints
                                   Write(6,'("Z=",i3," N=",i3," force=",a4," T=",f8.4,   &
                                &            " iq=",i2," l=",i2,                         &
                                &         " m=",i2," Qlm=",f10.5," -- RESTART i=",i6)')  &
                                           zNumber,nNumber,                              &
                                           forceVector(indexForce(indexJob)),            &
                                           temperatureVector(temperature),               &
                                           iq,lambdaJob(indexQjob(indexJob),iq),         &
                                           miuJob(indexQjob(indexJob),iq),               &
                                           qValueJob(indexQjob(indexJob),iq),            &
                                           indexJobVector_orig(iGrid,iCommon)
                                End Do
                             End If

                          End Do
                       End Do
                    End Do

                 End Do
              End Do
           End Do

           Deallocate(indexJobVector_orig,mapping_oldGrid,mapping_qlm,mapping_neck)

        Else

           iCommon = 0

           Allocate(mapping_oldGrid(1:numberQ,1:neck_numberPoints_orig))

           Do iz=1,proton_stepNumber
              Do in=1,neutron_stepNumber
                 Do force=1,counterForce
                    Do temperature=1,counterTemperature

                       iCommon=iCommon+1; iGrid=0

                       Do idxNeck=1,neck_numberPoints_orig
                          Do idxDefo=1,numberQ

                             iGrid=iGrid+1
                             mapping_oldGrid(idxDefo,idxNeck)=iGrid

                             indexJob_orig=indexJob_orig+1
                             indexJobVector_orig(iGrid,iCommon)=indexJob_orig

                          End Do
                       End Do

                    End Do
                 End Do
              End Do
           End Do

           iCommon=0; indexJob=0

           Do iz=1,proton_stepNumber
              zNumber = zVector(iz)
              If(x_proton_initValue.Gt.0.0D0) xzNumber=xzVector(iz)
              Do in=1,neutron_stepNumber
                 nNumber = nVector(in)
                 If(x_neutron_initValue.Gt.0.0D0) xnNumber=xnVector(in)
                 Do force=1,counterForce
                    Do temperature=1,counterTemperature

                       iCommon = iCommon + 1

                       Do i_neck=1,neck_numberPoints
                          Do ii=1,numberQ

                             indexJob=indexJob+1
                             indexQjob(indexJob)=ii

                             indexForce(indexJob)=force
                             indexTemperature(indexJob)=temperature
                             indexNeck(indexJob)=i_neck

                             IZ_JOB(indexJob)=zNumber
                             IN_JOB(indexJob)=nNumber

                             If(x_proton_initValue.Gt.0.0D0) Then
                                XIZJOB(indexJob)=xzNumber
                             Else
                                XIZJOB(indexJob)=Real(zNumber,Kind=Kind(1.0))
                             End If
                             If(x_neutron_initValue.Gt.0.0D0) Then
                                XINJOB(indexJob)=xnNumber
                             Else
                                XINJOB(indexJob)=Real(nNumber,Kind=Kind(1.0))
                             End If

                             idxDefo=mapping_qlm(ii)
                             idxNeck=mapping_neck(i_neck)
                             iGrid=mapping_oldGrid(idxDefo, idxNeck)

                             indexRestartJob(indexJob)=indexJobVector_orig(iGrid, iCommon)

                          End Do
                       End Do
                    End Do
                 End Do
              End Do
           End Do

           Deallocate(indexJobVector_orig,mapping_oldGrid,mapping_qlm,mapping_neck)

         End If

     End If

     If(modeMPI_deformation.Ge.3) Then

        NOFJOB = proton_stepNumber * neutron_stepNumber * counterForce * numberQ &
               * counterTemperature * neck_numberPoints

        Allocate(IZ_JOB(1:NOFJOB),IN_JOB(1:NOFJOB),indexQjob(1:NOFJOB),indexRestartJob(1:NOFJOB),&
                 indexForce(1:NOFJOB),indexTemperature(1:NOFJOB),indexNeck(1:NOFJOB))
        Allocate(XIZJOB(1:NOFJOB),XINJOB(1:NOFJOB)); XIZJOB(:)=-1.0D0; XINJOB(:)=-1.0D0

        nCommon = proton_stepNumber * neutron_stepNumber * counterForce * counterTemperature

        iCommon  = 0; indexJob = 0
        Do iz=1,proton_stepNumber
           zNumber = zVector(iz)
           If(x_proton_initValue.Gt.0.0D0) xzNumber=xzVector(iz)
           Do in=1,neutron_stepNumber
              nNumber = nVector(in)
              If(x_neutron_initValue.Gt.0.0D0) xnNumber=xnVector(in)
              Do force=1,counterForce

                 iCommon = iCommon + 1

                 Do temperature=1,counterTemperature
                    Do i_neck=1,neck_numberPoints
                       Do ii=1,numberQ

                          indexJob=indexJob+1
                          indexQjob(indexJob)=ii

                          indexForce(indexJob)=force
                          indexTemperature(indexJob)=temperature
                          indexNeck(indexJob)=i_neck

                          IZ_JOB(indexJob)=zNumber
                          IN_JOB(indexJob)=nNumber

                          If(x_proton_initValue.Gt.0.0D0) Then
                             XIZJOB(indexJob)=xzNumber
                          Else
                             XIZJOB(indexJob)=Real(zNumber,Kind=Kind(1.0))
                          End If
                          If(x_neutron_initValue.Gt.0.0D0) Then
                             XINJOB(indexJob)=xnNumber
                          Else
                             XINJOB(indexJob)=Real(nNumber,Kind=Kind(1.0))
                          End If

                          indexRestartJob(indexJob)=mapping_qlm_new(ii)

                          If(debug.Ge.4 .And. mpi_rank .Eq.0) Then
                             Do iq=1,numberConstraints
                                Write(6,'("Z=",i3," N=",i3," force=",a4," T=",f8.4,   &
                             &            " iq=",i2," l=",i2,                         &
                             &         " m=",i2," Qlm=",f10.5," -- RESTART i=",i6)')  &
                                        zNumber,nNumber,                              &
                                        forceVector(indexForce(indexJob)),            &
                                        temperatureVector(temperature),               &
                                        iq,lambdaJob(indexQjob(indexJob),iq),         &
                                        miuJob(indexQjob(indexJob),iq),               &
                                        qValueJob(indexQjob(indexJob),iq),            &
                                        ii
                             End Do
                          End If

                       End Do
                    End Do
                 End Do

              End Do
           End Do
        End Do

        Deallocate(mapping_qlm_new,mapping_neck)

     End If

     If(modeMPI_basis==1) Then

        Allocate(basis_size(1:basis_size_stepNumber))
        Allocate(basis_states(1:basis_states_stepNumber))
        Allocate(basis_freq(1:basis_frequency_stepNumber))
        Allocate(basis_defo(1:basis_deform_stepNumber))

        Do i_shel=1,basis_size_stepNumber
           basis_size(i_shel) = basis_size_initValue + (i_shel-1)*basis_size_stepValue
        End do
        Do i_stat=1,basis_states_stepNumber
           basis_states(i_stat) = basis_states_initValue + (i_stat-1)*basis_states_stepValue
        End do
        do i_freq=1,basis_frequency_stepNumber
           basis_freq(i_freq) = basis_frequency_initValue + Real(i_freq-1)*basis_frequency_stepValue
        End do
        do i_defo=1,basis_deform_stepNumber
           basis_defo(i_defo) = basis_deform_initValue + Real(i_defo-1)*basis_deform_stepValue
        End do

        NOFJOB = proton_stepNumber * neutron_stepNumber * numberQ   &
               * basis_frequency_stepNumber * basis_deform_stepNumber  &
               * basis_size_stepNumber * basis_states_stepNumber * counterForce

        Allocate(IZ_JOB(1:NOFJOB),IN_JOB(1:NOFJOB),indexQjob(1:NOFJOB),indexForce(1:NOFJOB),&
                 indexSizeJob(1:NOFJOB),indexStatesJob(1:NOFJOB),indexLengthJob(1:NOFJOB),  &
                 idxDeformJob(1:NOFJOB))
        Allocate(XIZJOB(1:NOFJOB),XINJOB(1:NOFJOB)); XIZJOB(:)=-1.0D0; XINJOB(:)=-1.0D0

        indexJob = 0

        Do iz=1,proton_stepNumber
           zNumber = zVector(iz)
           If(x_proton_initValue.Gt.0.0D0) xzNumber=xzVector(iz)
           Do in=1,neutron_stepNumber
              nNumber = nVector(in)
              If(x_neutron_initValue.Gt.0.0D0) xnNumber=xnVector(in)
              Do force=1,counterForce
                  Do ii=1,numberQ
                     Do i_shel=1,basis_size_stepNumber
                        Do i_stat=1,basis_states_stepNumber
                           Do i_freq=1,basis_frequency_stepNumber
                              Do i_defo=1,basis_deform_stepNumber

                                 indexJob = indexJob + 1

                                 IZ_JOB(indexJob)=zNumber
                                 IN_JOB(indexJob)=nNumber

                                 If(x_proton_initValue.Gt.0.0D0) Then
                                    XIZJOB(indexJob)=xzNumber
                                 Else
                                    XIZJOB(indexJob)=Real(zNumber,Kind=Kind(1.0))
                                 End If
                                 If(x_neutron_initValue.Gt.0.0D0) Then
                                    XINJOB(indexJob)=xnNumber
                                 Else
                                    XINJOB(indexJob)=Real(nNumber,Kind=Kind(1.0))
                                 End If

                                 indexQjob(indexJob)=ii
                                 indexForce(indexJob)=force

                                 indexSizeJob(indexJob)=i_shel
                                 indexStatesJob(indexJob)=i_stat
                                 indexLengthJob(indexJob)=i_freq
                                 idxDeformJob(indexJob)=i_defo

                              End Do
                           End Do
                        End Do
                     End Do
                  End Do
              End Do
           End Do
        End Do

     End If

  End Subroutine mpi_constructJobList

  !---------------------------------------------------------------------!
  !  Subroutine setting local (process-wise) HFODD variables based on   !
  !  existing list of tasks. Routine executed by each process.          !
  !---------------------------------------------------------------------!
  Subroutine mpi_getCurrentData(INDJOB,IN_FIX,IZ_FIX,SKYRME,GOGNAM, &
                                              INSIQN,IPSIQN,IDSIQN, &
                                              INSIQP,IPSIQP,IDSIQP, &
                                       INSIGN,IPSIGN,ISSIGN,IDSIGN, &
                                       INSIGP,IPSIGP,ISSIGP,IDSIGP, &
                                       ICONTI,IPCONT,ILCONT,IFCONT, &
                                                            I_GOGA, &
                                       SLOWEV,SLOWOD,SLOWPA,SLOWLI, &
                                                            IROTAT, &
                                                            NILXYZ, &
                                                            N_EXIT, &
                                                            NOITER, &
                                              FCHOM0,NOSCIL,NLIMIT, &
                                                     LIPKIN,LIPKIP, &
                                       NXHERM,NYHERM,NZHERM,IOPTGS, &
                                              IPAIRI,IPAHFB,IMFHFB, &
                                                            IWRILI, &
                                                            EPSITE, &
                                                            INIROT, &
                                                            IGAMMA, &
                                                            LANODD, &
                                       IFIBLN,INIBLN,IFIBLP,INIBLP, &
                                                            NMUCON, &
                                                            IFSHEL, &
                                                            IF_RPA, &
                                                            INDFIL, &
                                                            TEMP_T)

     Use hfodd_sizes

     Integer :: INDJOB,IN_FIX,IZ_FIX,INSIQN,IPSIQN,IDSIQN,INSIQP,IPSIQP,IDSIQP,&
                INSIGN,IPSIGN,ISSIGN,IDSIGN,INSIGP,IPSIGP,ISSIGP,IDSIGP,ICONTI,&
                IPCONT,ILCONT,IROTAT,NILXYZ,N_EXIT,NOITER,NOSCIL,NLIMIT,LIPKIN,&
                LIPKIP,NXHERM,NYHERM,NZHERM,IOPTGS,IPAIRI,IPAHFB,IMFHFB,IWRILI,&
                IFCONT,INIROT,IGAMMA,LANODD,IFIBLN,INIBLN,IFIBLP,INIBLP,NMUCON,&
                IFSHEL,IF_RPA,INDFIL,I_GOGA

     Integer :: IFLAGQ,IFLALQ,indexQ,ii,LAMBDA,MIU,IFNECK
     Integer :: i_shel,i_stat,i_freq,i_defo,force,temperature,i_neck,is_deformed
     Integer :: mpi_rank,mpi_err

     Character(Len=4) :: SKYRME,GOGNAM

     Real :: SLOWEV,SLOWOD,SLOWPA,SLOWLI,EPSITE,FCHOM0,STIFFQ, &
             QASKED,GALMUQ,QLINEA,Q0NECK,G_NECK,HBMASS,HBMRPA, &
             HBMINP,H_BARC,XMASSP,XMASSN,HBCOE2,HBAROX,HBAROY, &
             HBAROZ,HOMSCA,ALPHAR,TEMP_T,PIARGU,OVALUE,A_MASS, &
             HOMEGA,XINFIX,XIZFIX

     COMMON                                           &
            /FLOPAR/ XIZFIX,XINFIX
     COMMON                                           &
            /NCKVAL/ Q0NECK,G_NECK                    &
            /NCKFLA/ IFNECK
     COMMON                                           &
            /REALPH/ ALPHAR(0:NDLAMB,0:NDLAMB)
     COMMON                                           &
            /QCNSTR/ STIFFQ(0:NDMULT,-NDMULT:NDMULT), &
                     QASKED(0:NDMULT,-NDMULT:NDMULT), &
                     IFLAGQ(0:NDMULT,-NDMULT:NDMULT)
     COMMON                                           &
            /QLASTR/ GALMUQ(0:NDMULT,-NDMULT:NDMULT), &
                     QLINEA(0:NDMULT,-NDMULT:NDMULT), &
                     IFLALQ(0:NDMULT,-NDMULT:NDMULT)
     COMMON                                           &
            /PLANCK/ HBMASS,HBMRPA,HBMINP
     COMMON                                           &
            /PHYCON/ H_BARC,XMASSP,XMASSN,HBCOE2
     COMMON                                           &
            /BASPAR/ HBAROX,HBAROY,HBAROZ             &
            /SCALNG/ HOMSCA(NDKART,NDTWCE)

     Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err)

     ! Set proton, neutron number and Skyrme functional
     IN_FIX=IN_JOB(INDJOB); IZ_FIX=IZ_JOB(INDJOB)
     XINFIX=XINJOB(INDJOB); XIZFIX=XIZJOB(INDJOB)
     If(debug.Ge.1) Write(6,'("In mpi_getCurrentData - mpi_rank = ",i10," XN = ",f20.14," XZ = ",f20.14)') &
                                                       mpi_rank, XINFIX, XIZFIX
     If(I_GOGA.Gt.0) Then
        force=indexForce(INDJOB)
        GOGNAM=forceVector(force)
        SKYRME=GOGNAM
     Else
        force=indexForce(INDJOB)
        SKYRME=forceVector(force)
     End If

     ! Set temperature and index of the restart file (assumed to be in restart/ directory,
     ! and to have the generic form HFODD_XXXXXXXX.REC, where XXXXXXXX is a 8 digit integer
     If(modeMPI_deformation.Ge.1.Or.modeMPI_neck.Ge.1) Then
        temperature=indexTemperature(INDJOB)
        TEMP_T=temperatureVector(temperature)
        INDFIL=indexRestartJob(INDJOB)
        If(debug.Ge.1) Write(6,'("In mpi_getCurrentData - mpi_rank = ",i10," INDFIL = ",i10)') &
                                                          mpi_rank,INDFIL
     End If

     ! Set constraint on the neck, and active it
     If(modeMPI_neck.Ge.1) Then
        IFNECK=2
        i_neck=indexNeck(INDJOB)
        Q0NECK=neck_value(i_neck)
     End If

     ! Set constraints on multipole moments, and activate constraints
     indexQ=indexQjob(INDJOB)
     If(debug.Ge.1) Write(6,'("In mpi_getCurrentData - mpi_rank = ",i10," INDJOB = ",i10," indexQ = ",i10)') &
                                                       mpi_rank,INDJOB,indexQ
     Do ii=1,numberConstraints
        LAMBDA=Abs(lambdaJob(indexQ,ii))
        If(LAMBDA.Ge.NMUCON) NMUCON=LAMBDA
        MIU = miuJob(indexQ, ii)
        QASKED(LAMBDA,MIU) = qValueJob(indexQ, ii)
                                  IFLAGQ(LAMBDA,MIU) = 0
        If(turnConstraintOn.Eq.1) IFLAGQ(LAMBDA,MIU) = 1
                        IFLALQ(LAMBDA,MIU) = 0
        If(IF_RPA.Eq.1) IFLALQ(LAMBDA,MIU) =-1 ! RPA
        STIFFQ(LAMBDA,MIU) = 0.5D0
     End Do

     ! Set characteristics of the basis for this process
     If(modeMPI_basis.Eq.1) Then

        i_shel=indexSizeJob(INDJOB)
        i_stat=indexStatesJob(INDJOB)
        i_freq=indexLengthJob(INDJOB)
        i_defo=idxDeformJob(INDJOB)

        ! Number of shells and number of states
        NOSCIL = basis_size(i_shel)
        NLIMIT = basis_states(i_stat)

        ! Check if the requested configuration is deformed
        is_deformed = 0
        Do lambda=0,NDMULT
           Do miu=0,NDMULT
              If(QASKED(LAMBDA,MIU).Gt.1.D-12) Then
                 If(debug.Ge.4) Write(6,'("LAMBDA = ",i2," MIU = ",i2," QASKED = ",f20.14)') &
                                     LAMBDA,MIU,QASKED(LAMBDA,MIU)
                 is_deformed = 1
              End If
           End Do
        End Do
        ! Set frequency and deformation of the basis
        A_MASS=Real(IZ_FIX+IN_FIX); HOMEGA=41.0D0/A_MASS**(1.0D0/3.0D0)
        OVALUE=basis_freq(i_freq); FCHOM0=OVALUE/HOMEGA
        ALPHAR(2,0)=basis_defo(i_defo)
        If(debug.Ge.1) Write(6,'("NOSCIL=",i4," NLIMIT=",i4," FCHOM0=",f10.5," ALPHAR=",f10.5)') &
                             NOSCIL,NLIMIT,FCHOM0,ALPHAR(2,0)

     End If

   End Subroutine mpi_getCurrentData

   !---------------------------------------------------------------------!
   !  Subroutine that allows to redefine the requested values of the     !
   !  constraints based on the values of Qlm read from disk. This is     !
   !  used in restart mode only to effectively constrain previously free !
   !  collective variables.                                              !
   !  Example: Run 1 - Q20 grid, Q40 unconstrained                       !
   !           Run 2 - Same Q20 grid, new Q40 grid between 0 and 50.     !
   !                   Calling mpi_pathExploration will recenter the Q40 !
   !                   grid around the equilibrium values of Q40 from    !
   !                   Run 1.                                            !
   !---------------------------------------------------------------------!
   Subroutine mpi_pathExploration(IERROR,NMUCON)

     Use hfodd_sizes

     Integer :: IERROR,NMUCON,ITITLE,NFIPRI,IFLAGQ,LAMACT,MIUACT,i,lambda,miu

     Real :: QOLDIE,QSHIFT,STIFFQ,QASKED,QMUL_T,QMUL_P,QMUL_N

     COMMON                                           &
            /CFIPRI/ NFIPRI
     COMMON                                           &
            /QCNSTR/ STIFFQ(0:NDMULT,-NDMULT:NDMULT), &
                     QASKED(0:NDMULT,-NDMULT:NDMULT), &
                     IFLAGQ(0:NDMULT,-NDMULT:NDMULT)
     COMMON                                           &
            /QMULTI/ QMUL_N(0:NDMULT,-NDMULT:NDMULT), &
                     QMUL_P(0:NDMULT,-NDMULT:NDMULT), &
                     QMUL_T(0:NDMULT,-NDMULT:NDMULT)

     If(modeMPI_deformation.Eq.2.And.optimizePath.Eq.1.And.IERROR.Eq.0) Then

        ! Readjustement only for extra constraints not included in the original run
        Do i=1,numberConstraints
           LAMACT=qLambda_restart(i)
           MIUACT=qMiu_restart(i)
           QOLDIE=QASKED(LAMACT,MIUACT) ! Original value requested
           QSHIFT=QMUL_T(LAMACT,MIUACT) - qinitValue_restart(i) &
                                   -0.5D0*(qFinValue_restart(i) - qinitValue_restart(i))
           QASKED(LAMACT,MIUACT)=QSHIFT+QASKED(LAMACT,MIUACT) ! new value centered on unconstrained result
        End Do

        ! Displaying a message with the new values of the constraints
        Write(NFIPRI,'("*",77X,"*")')
        ITITLE=0
        Do lambda=1,NMUCON
           Do miu=-lambda,lambda
              If(IFLAGQ(lambda,miu).Eq.1) Then
                 If(ITITLE.Eq.0) Then
                     Write(NFIPRI,'(79("*"),/,"*",77X,"*")')
                     Write(NFIPRI,'("*",1X,"UPDATED MULTIPOLE CONSTRAINTS: ",&
                    &     "LAMBDA=",I2,"  MIU=",I2,2X,"MOMENT=",F7.3,12X,"*")') &
                           lambda,miu,QASKED(lambda,miu)
                     ITITLE=1
                 Else
                     Write(NFIPRI,'("*",32X,"LAMBDA=",I2,"  MIU=",I2,2X,"MOMENT=",F7.3,12X,"*")') &
                           lambda,miu,QASKED(lambda,miu)
                 End If
              End If
           End Do
        End Do

     End If

   End Subroutine mpi_pathExploration

   !---------------------------------------------------------------------!
   !  Subroutine reading the HFODD record file. It is a modified version !
   !  of routine RECORD(), where storage of big arrays is made through   !
   !  allocatable Fortran arrays.                                        !
   !---------------------------------------------------------------------!
   Subroutine read_hfodd(NFIREP,FILREP,NUMITE,NXHERM,NYHERM,NZHERM,IPCONT, &
                         ILCONT,IACONT,IMCONT,IRCONT,IRENMA,IRENIN,LIPKIN, &
                         LIPKIP,REFERN,REFERP,REDELN,REDELP,REFE2N,REFE2P, &
                         IFIBLN,INIBLN,IFIBLP,INIBLP,NMUMAX,NSIMAX,NMUCON, &
                                                            ISHIFT,IERROR)

     Use hfodd_sizes

     Integer :: currentVersion, minimalVersion, mpi_rank, mpi_err, iq
     Integer :: IX,IY,IZ,K,L,LDBASE,NFIPRI,IFNECK,IFIRED,LAMBDA,MIU, &
                NFIREP,NUMITE,NXHERM,NYHERM,NZHERM,MDMULT,IPCONT,    &
                ILCONT,IACONT,IMCONT,IRCONT,IRENMA,IRENIN,LIPKIN,    &
                LIPKIP,IFIBLN,INIBLN,IFIBLP,INIBLP,NMUMAX,NSIMAX,    &
                NMUCON,ISHIFT,IERROR,IVERIN,MXHERM,MYHERM,MZHERM,    &
                MMUMAX,MSIMAX,MDBASE,ICHARG,IREVER,KARTEZ,IBASE,ii
     Integer, Allocatable :: IFLALQ(:,:),IFLALS(:,:),IDEBLO(:)

     character(Len=68) :: FILREP

     Real :: REFERN,REFERP,REDELN,REDELP,REFE2N,REFE2P,Q0NECK,G_NECK,DUMMY
     Real, Allocatable :: ANGU_N(:),ANGU_P(:),ANGU_T(:),SPIN_N(:),SPIN_P(:),SPIN_T(:)
     Real, Allocatable :: QMUL_N(:,:),QMUL_P(:,:),QMUL_T(:,:),SMUL_N(:,:),SMUL_P(:,:), &
                          SMUL_T(:,:),GALMUQ(:,:),QLINEA(:,:),GALMUS(:,:),SLINEA(:,:), &
                          RALMUQ(:,:),RALMUS(:,:)
     Real, Allocatable :: FERINI(:),DELINI(:),FE2INI(:),RSHIFT(:),HBMREN(:),ROTREN(:), &
                                      HBMREA(:),ROTREA(:)
     Real, Allocatable :: DLINSN(:),DLINSP(:),DLINST(:),ELINSN(:),ELINSP(:),ELINST(:), &
                          TLINSN(:),TLINSP(:),TLINST(:),ALINLN(:),ALINLP(:),ALINLT(:), &
                          DROTSN(:),DROTSP(:),DROTST(:),EROTSN(:),EROTSP(:),EROTST(:), &
                          TROTSN(:),TROTSP(:),TROTST(:),AROTLN(:),AROTLP(:),AROTLT(:)
     Real, Allocatable :: VN_MAS(:,:,:),VN_CEN(:,:,:),VP_MAS(:,:,:),VP_CEN(:,:,:), &
                          DPRRHO(:,:,:),DNRRHO(:,:,:)
     Real, Allocatable :: VN_KIS(:,:,:,:),VN_SPI(:,:,:,:),VN_CUR(:,:,:,:), &
                          VP_KIS(:,:,:,:),VP_SPI(:,:,:,:),VP_CUR(:,:,:,:)
     Real, Allocatable :: VN_SOR(:,:,:,:,:),VP_SOR(:,:,:,:,:)

     Complex :: C_ZERO
     Complex, Allocatable :: DN_RHO(:,:,:),DN_TAU(:,:,:),DN_LPR(:,:,:),DN_DIV(:,:,:), &
                             DP_RHO(:,:,:),DP_TAU(:,:,:),DP_LPR(:,:,:),DP_DIV(:,:,:),&
                             WN_CEN(:,:,:),WP_CEN(:,:,:),WAVBLO(:,:,:)
     Complex, Allocatable :: DN_SPI(:,:,:,:),DN_KIS(:,:,:,:),DN_GRR(:,:,:,:),DN_LPS(:,:,:,:),&
                             DN_ROS(:,:,:,:),DN_ROC(:,:,:,:),DN_CUR(:,:,:,:),                &
                             DP_SPI(:,:,:,:),DP_KIS(:,:,:,:),DP_GRR(:,:,:,:),DP_LPS(:,:,:,:),&
                             DP_ROS(:,:,:,:),DP_ROC(:,:,:,:),DP_CUR(:,:,:,:)
     Complex, Allocatable :: DN_SCU(:,:,:,:,:),DN_DES(:,:,:,:,:),DP_SCU(:,:,:,:,:),DP_DES(:,:,:,:,:)

     COMMON                 &
            /DIMENS/ LDBASE
     COMMON                 &
            /CFIPRI/ NFIPRI
     COMMON                 &
            /NCKFLA/ IFNECK

     C_ZERO=Cmplx(0.0D0,0.0D0)

     currentVersion = 16; minimalVersion =  4

     Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err)

     Open(UNIT=NFIREP,FILE=FILREP,STATUS='OLD',FORM='UNFORMATTED',IOSTAT=IERROR)
     If(IERROR.Ne.0) Then
         Write(NFIPRI,'("read_hfodd() - Error in opening file ",i5)') NFIREP
         Return
     End If

     IFIRED=0

     Read(NFIREP,IOSTAT=IERROR) IVERIN
     If(IERROR.Ne.0) Then
         Write(NFIPRI,'("read_hfodd() - error in reading IVERIN")')
         Return
     End If

     If(.Not.Allocated(WN_CEN)) &
        Allocate(WN_CEN(NDXHRM,NDYHRM,NDZHRM),VN_MAS(NDXHRM,NDYHRM,NDZHRM),VN_CEN(NDXHRM,NDYHRM,NDZHRM),&
                 WP_CEN(NDXHRM,NDYHRM,NDZHRM),VP_MAS(NDXHRM,NDYHRM,NDZHRM),VP_CEN(NDXHRM,NDYHRM,NDZHRM))
     If(.Not.Allocated(VN_KIS)) &
        Allocate(VN_KIS(NDXHRM,NDYHRM,NDZHRM,NDKART),VN_SPI(NDXHRM,NDYHRM,NDZHRM,NDKART),&
                 VN_CUR(NDXHRM,NDYHRM,NDZHRM,NDKART),VP_KIS(NDXHRM,NDYHRM,NDZHRM,NDKART),&
                 VP_SPI(NDXHRM,NDYHRM,NDZHRM,NDKART),VP_CUR(NDXHRM,NDYHRM,NDZHRM,NDKART))
     If(.Not.Allocated(VN_SOR)) &
        Allocate(VN_SOR(NDXHRM,NDYHRM,NDZHRM,NDKART,NDKART),VP_SOR(NDXHRM,NDYHRM,NDZHRM,NDKART,NDKART))

     If(.Not.Allocated(DN_RHO)) &
        Allocate(DN_RHO(NDXHRM,NDYHRM,NDZHRM),DN_TAU(NDXHRM,NDYHRM,NDZHRM),DN_LPR(NDXHRM,NDYHRM,NDZHRM),&
                 DN_DIV(NDXHRM,NDYHRM,NDZHRM),DP_RHO(NDXHRM,NDYHRM,NDZHRM),DP_TAU(NDXHRM,NDYHRM,NDZHRM),&
                 DP_LPR(NDXHRM,NDYHRM,NDZHRM),DP_DIV(NDXHRM,NDYHRM,NDZHRM))
     If(.Not.Allocated(DN_SPI)) &
        Allocate(DN_SPI(NDXHRM,NDYHRM,NDZHRM,NDKART),DN_KIS(NDXHRM,NDYHRM,NDZHRM,NDKART),&
                 DN_GRR(NDXHRM,NDYHRM,NDZHRM,NDKART),DN_LPS(NDXHRM,NDYHRM,NDZHRM,NDKART),&
                 DN_ROS(NDXHRM,NDYHRM,NDZHRM,NDKART),DN_ROC(NDXHRM,NDYHRM,NDZHRM,NDKART),&
                 DN_CUR(NDXHRM,NDYHRM,NDZHRM,NDKART))
     If(.Not.Allocated(DP_SPI)) &
        Allocate(DP_SPI(NDXHRM,NDYHRM,NDZHRM,NDKART),DP_KIS(NDXHRM,NDYHRM,NDZHRM,NDKART),&
                 DP_GRR(NDXHRM,NDYHRM,NDZHRM,NDKART),DP_LPS(NDXHRM,NDYHRM,NDZHRM,NDKART),&
                 DP_ROS(NDXHRM,NDYHRM,NDZHRM,NDKART),DP_ROC(NDXHRM,NDYHRM,NDZHRM,NDKART),&
                 DP_CUR(NDXHRM,NDYHRM,NDZHRM,NDKART))
     If(.Not.Allocated(DN_SCU)) &
        Allocate(DN_SCU(NDXHRM,NDYHRM,NDZHRM,NDKART,NDKART),DN_DES(NDXHRM,NDYHRM,NDZHRM,NDKART,NDKART),&
                 DP_SCU(NDXHRM,NDYHRM,NDZHRM,NDKART,NDKART),DP_DES(NDXHRM,NDYHRM,NDZHRM,NDKART,NDKART))

     If(.Not.Allocated(ANGU_N)) &
        Allocate(ANGU_N(1:NDKART),ANGU_P(1:NDKART),ANGU_T(1:NDKART),&
                 SPIN_N(1:NDKART),SPIN_P(1:NDKART),SPIN_T(1:NDKART))

     If(.Not.Allocated(QMUL_N)) &
        Allocate(QMUL_N(0:NDMULT,-NDMULT:NDMULT),QMUL_P(0:NDMULT,-NDMULT:NDMULT),QMUL_T(0:NDMULT,-NDMULT:NDMULT),&
                 SMUL_N(0:NDMULT,-NDMULT:NDMULT),SMUL_P(0:NDMULT,-NDMULT:NDMULT),SMUL_T(0:NDMULT,-NDMULT:NDMULT),&
                 GALMUQ(0:NDMULT,-NDMULT:NDMULT),QLINEA(0:NDMULT,-NDMULT:NDMULT),RALMUQ(0:NDMULT,-NDMULT:NDMULT),&
                 GALMUS(0:NDMULT,-NDMULT:NDMULT),SLINEA(0:NDMULT,-NDMULT:NDMULT),RALMUS(0:NDMULT,-NDMULT:NDMULT))
     If(.Not.Allocated(IFLALQ)) &
        Allocate(IFLALQ(0:NDMULT,-NDMULT:NDMULT),IFLALS(0:NDMULT,-NDMULT:NDMULT))

     If(.Not.Allocated(FERINI)) &
        Allocate(FERINI(0:NDISOS),DELINI(0:NDISOS),FE2INI(0:NDISOS),RSHIFT(1:NDKART),&
                 HBMREN(1:NDKART),ROTREN(1:NDKART),HBMREA(1:NDKART),ROTREA(1:NDKART))

     If(.Not.Allocated(DLINSN)) &
        Allocate(DLINSN(0:NDKART),DLINSP(0:NDKART),DLINST(0:NDKART), &
                 ELINSN(0:NDKART),ELINSP(0:NDKART),ELINST(0:NDKART), &
                 TLINSN(0:NDKART),TLINSP(0:NDKART),TLINST(0:NDKART), &
                 ALINLN(0:NDKART),ALINLP(0:NDKART),ALINLT(0:NDKART), &
                 DROTSN(0:NDKART),DROTSP(0:NDKART),DROTST(0:NDKART), &
                 EROTSN(0:NDKART),EROTSP(0:NDKART),EROTST(0:NDKART), &
                 TROTSN(0:NDKART),TROTSP(0:NDKART),TROTST(0:NDKART), &
                 AROTLN(0:NDKART),AROTLP(0:NDKART),AROTLT(0:NDKART))

     If(.Not.Allocated(WAVBLO)) &
        Allocate(WAVBLO(1:2*NDBASE,0:NDISOS,0:NDREVE))
     If(.Not.Allocated(IDEBLO)) &
        Allocate(IDEBLO(0:NDISOS))

     If(IVERIN.Ge.minimalVersion.AND.IVERIN.Le.currentVersion) Then

        Read(NFIREP,IOSTAT=IERROR) NUMITE
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading NUMITE")')
           Return
        End If

        Read(NFIREP,IOSTAT=IERROR) MXHERM,MYHERM,MZHERM
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading MXHERM")')
           Return
        End If

        If(MXHERM.Ne.NXHERM.OR.MYHERM.Ne.NYHERM.OR.MZHERM.Ne.NZHERM) Then
           Write(NFIPRI,'(1X,19("/"),"  WRONG NXHERM,NYHERM,NZHERM:",3I3,2X,19("/"))') MXHERM,MYHERM,MZHERM
           Stop '  WRONG NXHERM,NYHERM,NZHERM IN RECORD'
        End If

        Read(NFIREP,IOSTAT=IERROR) (((VN_MAS(IX,IY,IZ),VN_CEN(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM)
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading VN_MAS")')
           Return
        End If

        Read(NFIREP,IOSTAT=IERROR) ((((VN_KIS(IX,IY,IZ,L),VN_SPI(IX,IY,IZ,L),VN_CUR(IX,IY,IZ,L), &
                                               IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM),L=1,NDKART)
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading VN_KIS")')
           Return
        End If

        Read(NFIREP,IOSTAT=IERROR) (((((VN_SOR(IX,IY,IZ,L,K),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM),L=1,NDKART),&
                                                                                                    K=1,NDKART)
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading VN_SOR")')
           Return
        End If

        Read(NFIREP,IOSTAT=IERROR) (((VP_MAS(IX,IY,IZ),VP_CEN(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM)
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading VP_MAS")')
           Return
        End If

        Read(NFIREP,IOSTAT=IERROR) ((((VP_KIS(IX,IY,IZ,L),VP_SPI(IX,IY,IZ,L),VP_CUR(IX,IY,IZ,L), &
                                               IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM),L=1,NDKART)
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading VP_KIS")')
           Return
        End If

        Read(NFIREP,IOSTAT=IERROR) (((((VP_SOR(IX,IY,IZ,L,K),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM),L=1,NDKART),&
                                                                                                    K=1,NDKART)
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading VP_SOR")')
           Return
        End If

        If(IVERIN.Le.6) Then

           If(IVERIN.Gt.5) Then
               Read(NFIREP,IOSTAT=IERROR) (((DNRRHO(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM)
               If(IERROR.Ne.0) Then
                   Write(NFIPRI,'("read_hfodd() - Error in reading DNRRHO")')
                   Return
               End If
           Else
               DNRRHO(:,:,:)=0.0D0
           End If

           Read(NFIREP,IOSTAT=IERROR) (((DPRRHO(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM)
           If(IERROR.Ne.0) Then
              Write(NFIPRI,'("read_hfodd() - Error in reading DPRRHO")')
              Return
           End If
           DN_RHO(:,:,:)=DNRRHO(:,:,:)*cmplx(1.0d0,0.0d0)
           DP_RHO(:,:,:)=DPRRHO(:,:,:)*cmplx(1.0d0,0.0d0)

        Else

           Read(NFIREP,IOSTAT=IERROR) (((DN_RHO(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM)
           If(IERROR.Ne.0) Then
              Write(NFIPRI,'("read_hfodd() - Error in reading DN_RHO")')
              Return
           End If

           Read(NFIREP,IOSTAT=IERROR)(((DP_RHO(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM)
           If(IERROR.Ne.0) Then
              Write(NFIPRI,'("read_hfodd() - Error in reading DP_RHO")')
              Return
           End If

        End If

        Read(NFIREP,IOSTAT=IERROR) MMUMAX
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading NMUMAX")')
           Return
        End If

        If(MMUMAX.Lt.NMUMAX) Then

           Write(NFIPRI,'(1X,20("/"),"  WARNING:  TOO FEW MULTIPOLE MOMENTS",1X,20("/"))')
           Write(NFIPRI,'(1X,20("/"),"  FOUND ON THE REPLAY FILE.          ",1X,20("/"))')
           Write(NFIPRI,'(1X,20("/"),"  REQUESTED MAXIMUM  MULTIPOLARITY=",I2,1X,20("/"))') NMUMAX
           Write(NFIPRI,'(1X,20("/"),"  FOUND     MAXIMUM  MULTIPOLARITY=",I2,1X,20("/"))') MMUMAX
           Write(NFIPRI,'(1X,20("/"),"  A CONSTRAINT ON A MISSING MULTIPOLE",1X,20("/"))')
           Write(NFIPRI,'(1X,20("/"),"  WILL BE  INCLUDED  ONLY  AFTER  THE",1X,20("/"))')
           Write(NFIPRI,'(1X,20("/"),"  FIRST DIAGONALIZATION IS  PERFORMED",1X,20("/"))')

           Do LAMBDA=MMUMAX+1,NMUMAX
              Do MIU=-LAMBDA,LAMBDA
                 QMUL_N(LAMBDA,MIU)=0.0D0
                 QMUL_P(LAMBDA,MIU)=0.0D0
                 QMUL_T(LAMBDA,MIU)=0.0D0
              End Do
           End Do

        End If

        Read(NFIREP,IOSTAT=IERROR) ((QMUL_N(LAMBDA,MIU),QMUL_P(LAMBDA,MIU),QMUL_T(LAMBDA,MIU), &
                                                    MIU=-LAMBDA,LAMBDA),LAMBDA=0,MMUMAX)
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading QMUL_N")')
           Return
        End If

        Do LAMBDA=0,NSIMAX
           Do MIU=-LAMBDA,LAMBDA
              SMUL_N(LAMBDA,MIU)=0.0D0
              SMUL_P(LAMBDA,MIU)=0.0D0
              SMUL_T(LAMBDA,MIU)=0.0D0
           End Do
        End Do

        If(IVERIN.Gt.4) Then

           Read(NFIREP,IOSTAT=IERROR) MSIMAX
           If(IERROR.Ne.0) Then
              Write(NFIPRI,'("read_hfodd() - Error in reading MSIMAX")')
              Return
           End If

           If(MSIMAX.Lt.NSIMAX) Then

              Write(NFIPRI,'(1X,20(1H/),"  WARNING:  TOO FEW SURFACE   MOMENTS",1X,20(1H/))')
              Write(NFIPRI,'(1X,20(1H/),"  FOUND ON THE REPLAY FILE.          ",1X,20(1H/))')
              Write(NFIPRI,'(1X,20(1H/),"  REQUESTED MAXIMUM  MULTIPOLARITY=",I2,1X,20(1H/))') NSIMAX
              Write(NFIPRI,'(1X,20(1H/),"  FOUND     MAXIMUM  MULTIPOLARITY=",I2,1X,20(1H/))') MSIMAX
              Write(NFIPRI,'(1X,20(1H/),"  A CONSTRAINT ON A MISSING MULTIPOLE",1X,20(1H/))')
              Write(NFIPRI,'(1X,20(1H/),"  WILL BE  INCLUDED  ONLY  AFTER  THE",1X,20(1H/))')
              Write(NFIPRI,'(1X,20(1H/),"  FIRST DIAGONALIZATION IS  PERFORMED",1X,20(1H/))')

              Do LAMBDA=MSIMAX+1,NSIMAX
                 Do MIU=-LAMBDA,LAMBDA
                    SMUL_N(LAMBDA,MIU)=0.0D0
                    SMUL_P(LAMBDA,MIU)=0.0D0
                    SMUL_T(LAMBDA,MIU)=0.0D0
                 End Do
              End Do

           End If

           Read(NFIREP,IOSTAT=IERROR) ((SMUL_N(LAMBDA,MIU),SMUL_P(LAMBDA,MIU),SMUL_T(LAMBDA,MIU), &
                                              MIU=-LAMBDA,LAMBDA),LAMBDA=0,MSIMAX)
           If(IERROR.Ne.0) Then
              Write(NFIPRI,'("read_hfodd() - Error in reading SMUL_N")')
              Return
           End If

        End If

        Read(NFIREP,IOSTAT=IERROR) ISHIFT,(RSHIFT(L),L=1,NDKART)
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading ISHIFT")')
           Return
        End If

        Read(NFIREP,IOSTAT=IERROR) (ANGU_N(L),ANGU_P(L),ANGU_T(L),SPIN_N(L),SPIN_P(L),SPIN_T(L),L=1,NDKART)
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading ANGU_N")')
           Return
        End If

        Read(NFIREP,IOSTAT=IERROR) (((WN_CEN(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM)
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading WN_CEN")')
           Return
        End If

        Read(NFIREP,IOSTAT=IERROR) (((WP_CEN(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM)
        If(IERROR.Ne.0) Then
          Write(NFIPRI,'("read_hfodd() - Error in reading WP_CEN")')
          Return
        End If

        Read(NFIREP,IOSTAT=IERROR) REFERN,REDELN,REFERP,REDELP
        If(IERROR.Ne.0) Then
           Write(NFIPRI,'("read_hfodd() - Error in reading REFERN")')
           Return
        End If

        ! LN lambda_2 (IVERIN=8)
        If(IVERIN.Gt.7) Then
           Read(NFIREP,IOSTAT=IERROR) REFE2N,REFE2P
           If(IERROR.Ne.0) Then
              Write(NFIPRI,'("read_hfodd() - Error in reading REFE2N")')
              Return
           End If
        Else
           REFE2N=FE2INI(0); REFE2P=FE2INI(1)
           If((LIPKIN.EQ.1.OR.LIPKIP.EQ.1).AND.ILCONT.EQ.1) Then
               Write(NFIPRI,'(1X,13("/")," WARNING: THE LIPKIN-NOGAMI LAMBDA2S HAVE NOT BEEN",1X,13("/"))')
               Write(NFIPRI,'(1X,13("/")," FOUND ON THE  OLD  REPLAY  FILE  VERSION  NO. =",I2,1X,13("/"))') IVERIN
               Write(NFIPRI,'(1X,13("/")," THE  LN  CORRECTIONS CANNOT BE SMOOTHLY CONTINUED",1X,13("/"))')
           End If
        End If

        ! Blocked state wave-function (IVERIN=9)
        IDEBLO(0)=0; IDEBLO(1)=0
        If(IVERIN.Gt.8) Then

           Read(NFIREP,IOSTAT=IERROR) MDBASE,(IDEBLO(ICHARG),ICHARG=0,NDISOS)
           If(IERROR.Ne.0) Then
              Write(NFIPRI,'("read_hfodd() - Error in reading MDBASE")')
              Return
           End If

           If(MDBASE.Le.NDBASE) Then
              Read(NFIREP,IOSTAT=IERROR) (((WAVBLO(IBASE,ICHARG,IREVER),IBASE=1,MDBASE),ICHARG=0,NDISOS),&
                                                                IREVER=0,NDREVE)
              If(IERROR.Ne.0) Then
                 Write(NFIPRI,'("read_hfodd() - Error in reading WAVBLO")')
                 Return
              End If
           Else
              Read(NFIREP,IOSTAT=IERROR)
              If(IERROR.Ne.0) Then
                 Write(NFIPRI,'("read_hfodd() - Error in reading ...nothing...")')
                 Return
              End If
           End If

           If(((IFIBLN.EQ.1.AND.INIBLN.EQ.0).OR.(IFIBLP.EQ.1.AND.INIBLP.EQ.0)).AND. &
                (MDBASE.Ne.LDBASE.OR.MDBASE.Gt.NDBASE)) Then
                 Write(NFIPRI,'(1X,10("/")," INCOMPATIBLE DIMENSIONS FOUND ON THE RECORD FILE:",1X,10("/"))')
                 Write(NFIPRI,'(1X,10("/")," LDBASE=",I4," MDBASE=",I4," NDBASE=",I4,1X,10("/"))') LDBASE,MDBASE,NDBASE
                 Stop 'MDBASE<>LDBASE OR MDBASE>NDBASE IN RECORD'
           End If

        End If

        If(IVERIN.Ge.10) Then
           Read(NFIREP,IOSTAT=IERROR) (ALINLT(KARTEZ),KARTEZ=0,NDKART)
           If(IERROR.Ne.0) Then
              Write(NFIPRI,'("read_hfodd() - Error in reading ALINLT")')
              Return
           End If
        Else
           ALINLT(:)=0.0D0
           If(IMCONT.EQ.1) Then
              Write(NFIPRI,'(1X,13(1H/)," WARNING: THE AVERAGE LINEAR MOMENTA HAVE NOT BEEN",1X,13(1H/))')
              Write(NFIPRI,'(1X,13(1H/)," FOUND ON THE  OLD  REPLAY  FILE  VERSION  NO. =",I2,1X,13(1H/))') IVERIN
              Write(NFIPRI,'(1X,13(1H/)," THE EXACT CM CORREC. CANNOT BE SMOOTHLY CONTINUED",1X,13(1H/))')
           End If
        End If

        If(IVERIN.Ge.11) Then
           Read(NFIREP,IOSTAT=IERROR) (AROTLT(KARTEZ),KARTEZ=0,NDKART)
           If(IERROR.Ne.0) Then
              Write(NFIPRI,'("read_hfodd() - Error in reading AROTLT")')
              Return
           End If
        Else
           AROTLT(:)=0.0D0
           If(IRCONT.Eq.1) Then
              Write(NFIPRI,'(1X,13(1H/)," WARNING: THE AVERAGE ANGULAR MOMENTA HAVE NOT BEEN",1X,13(1H/))')
              Write(NFIPRI,'(1X,13(1H/)," FOUND ON  THE  OLD  REPLAY  FILE  VERSION  NO. =",I2,1X,13(1H/))') IVERIN
              Write(NFIPRI,'(1X,13(1H/)," THE EXACT ROT CORREC. CANNOT BE SMOOTHLY CONTINUED",1X,13(1H/))')
           End If
        End If

        If(IVERIN.Ge.12) Then

           Read(NFIREP,IOSTAT=IERROR) MDMULT
           If(IERROR.Ne.0) Then
              Write(NFIPRI,'("read_hfodd() - Error in reading MDMULT")')
              Return
           End If

           RALMUQ(:,:)=0.0D0; RALMUS(:,:)=0.0D0

           If(MDMULT.Ne.NDMULT) Then

              Write(NFIPRI,'(1X,20(1H/),"  WARNING: NOT ALL LINEAR CONSTRAINTS",1X,20(1H/))')
              Write(NFIPRI,'(1X,20(1H/),"  FOUND ON THE REPLAY FILE.          ",1X,20(1H/))')
              Write(NFIPRI,'(1X,20(1H/),"  REQUESTED MAXIMUM  MULTIPOLARITY=",I2,1X,20(1H/))') NDMULT
              Write(NFIPRI,'(1X,20(1H/),"  FOUND     MAXIMUM  MULTIPOLARITY=",I2,1X,20(1H/))') MDMULT
              Write(NFIPRI,'(1X,20(1H/),"  A CONSTRAINT ON A MISSING MULTIPOLE",1X,20(1H/))')
              Write(NFIPRI,'(1X,20(1H/),"  WILL BE  INCLUDED  ONLY  AFTER  THE",1X,20(1H/))')
              Write(NFIPRI,'(1X,20(1H/),"  FIRST DIAGONALIZATION IS  PERFORMED",1X,20(1H/))')

              RALMUQ(:,:)=QLINEA(:,:); RALMUS(:,:)=SLINEA(:,:)

              Read(NFIREP,IOSTAT=IERROR) ((DUMMY,DUMMY,MIU=-LAMBDA,LAMBDA),LAMBDA=0,MDMULT)
              If(IERROR.Ne.0) Then
                 Write(NFIPRI,'("read_hfodd() - Error in reading DUMMY")')
                 Return
              End If

           Else

               Read(NFIREP,IOSTAT=IERROR) ((RALMUQ(LAMBDA,MIU),RALMUS(LAMBDA,MIU),MIU=-LAMBDA,LAMBDA),&
                                                               LAMBDA=0,MDMULT)
               If(IERROR.Ne.0) Then
                  Write(NFIPRI,'("read_hfodd() - Error in reading RALMUQ")')
                  Return
               End If

               ! If not allocated, allocate; if already allocated, deallocate first and re-allocate (
               ! (enables vayring number of constraints)
               If(Allocated(flex_lambda)) Deallocate(flex_lambda,flex_miu,flex_qlm)
               Allocate(flex_lambda(numberConstraints),flex_miu(numberConstraints),flex_qlm(numberConstraints))
               Do iq=1,numberConstraints
                  lambda = qLambda(iq); miu = qMiu(iq)
                  flex_lambda(iq)=lambda
                  flex_miu(iq)=miu
                  flex_qlm(iq)=QMUL_T(lambda,miu)
               End Do

           End If

        Else

            RALMUQ(:,:)=QLINEA(:,:); RALMUS(:,:)=SLINEA(:,:)

            If(IACONT.EQ.1) Then
               Write(NFIPRI,'(1X,13(1H/)," WARNING: THE LINEAR CONSTRAINTS  H A V E  NOT BEEN",1X,13(1H/))')
               Write(NFIPRI,'(1X,13(1H/)," FOUND ON  THE  OLD  REPLAY  FILE  VERSION  NO. =",I2,1X,13(1H/))') IVERIN
               Write(NFIPRI,'(1X,13(1H/)," THE EXACT AUGM.LAGRA. CANNOT BE SMOOTHLY CONTINUED",1X,13(1H/))')
            End If

        End If

        If(IVERIN.Ge.13) Then

            Read(NFIREP,IOSTAT=IERROR) (HBMREA(KARTEZ),KARTEZ=1,NDKART)
            If(IERROR.Ne.0) Then
               Write(NFIPRI,'("read_hfodd() - Error in reading HBMREA")')
               Return
            End If

            Read(NFIREP,IOSTAT=IERROR) (ROTREA(KARTEZ),KARTEZ=1,NDKART)
            If(IERROR.Ne.0) Then
               Write(NFIPRI,'("read_hfodd() - Error in reading ROTREA")')
               Return
            End If

            If(IMCONT.EQ.1.AND.IRENMA.Ge.1) Then
               HBMREN(:)=HBMREA(:)
            End If

            If(IRCONT.EQ.1.AND.IRENIN.Ge.1) Then
               ROTREN(:)=ROTREA(:)
            End If

        Else

            If(IMCONT.EQ.1.AND.IRENMA.Ge.1) Then
               Write(NFIPRI,'(1X,13(1H/)," WARNING:  THE  RENORMALIZED MASSES  HAVE NOT BEEN",1X,13(1H/))')
               Write(NFIPRI,'(1X,13(1H/)," FOUND ON  THE  OLD  REPLAY FILE  VERSION  NO. =",I2,1X,13(1H/))') IVERIN
               Write(NFIPRI,'(1X,13(1H/)," THE EXACT CM CORREC. CANNOT BE SMOOTHLY CONTINUED",1X,13(1H/))')
            End If

            If(IRCONT.EQ.1.AND.IRENIN.Ge.1) Then
               Write(NFIPRI,'(1X,13(1H/)," WARNING:  THE   RENORMALIZED INERTIA HAVE NOT BEEN",1X,13(1H/))')
               Write(NFIPRI,'(1X,13(1H/)," FOUND  ON  THE  OLD  REPLAY FILE  VERSION  NO. =",I2,1X,13(1H/))') IVERIN
               Write(NFIPRI,'(1X,13(1H/)," THE EXACT ROT CORREC. CANNOT BE SMOOTHLY CONTINUED",1X,13(1H/))')
            End If

        End If

        ! Linear Constraints for the neck (IVERIN=14)
        If(IVERIN.Ge.14) Then
            Read(NFIREP,IOSTAT=IERROR) G_NECK
            If(IERROR.Ne.0) Then
               Write(NFIPRI,'("read_hfodd() - Error in reading MDMULT")')
               Return
            End If
        Else
            G_NECK=0.0D0
            If(IFNECK.Ge.2) Then
               Write(NFIPRI,'(1X,13(1H/)," WARNING : THE  NECK  CONSTRAINTS  H A S  NOT  BEEN",1X,13(1H/))')
               Write(NFIPRI,'(1X,13(1H/)," FOUND ON  THE  OLD  REPLAY  FILE  VERSION  NO. =",I2,1X,13(1H/))') IVERIN
               Write(NFIPRI,'(1X,13(1H/)," THE NECK CONSTRAINT CAN NOT BE  SMOOTHLY CONTINUED",1X,13(1H/))')
            End If
        End If

        IFIRED=1

     End If

     If(IFIRED.EQ.0) Then
        Write(NFIPRI,'(/,1X,10(1H/),"  CANNOT read THE REPLAY FILE VERSION:",I3,2X,10(1H/),/)') IVERIN
        Stop '  WRONG IVERIN IN RECORD'
     End If

     If(IFIBLN.EQ.1.AND.INIBLN.EQ.0.AND.IDEBLO(0).Ne.1) Then
        Write(NFIPRI,'(1X,10(1H/)," NEUTRON SINGLE-PARTICLE WF NOT FOUND ON THE RECORD FILE",1X,10(1H/))')
        Write(NFIPRI,'(1X,10(1H/)," RERUN THE BLOCKING CALCULATION WITH INIBLN=1           ",1X,10(1H/))')
        Stop ' NEUTRON SINGLE-PARTICLE WF NOT FOUND IN RECORD'
     End If

     If(IFIBLP.EQ.1.AND.INIBLP.EQ.0.AND.IDEBLO(1).Ne.1) Then
        Write(NFIPRI,'(1X,10(1H/)," PROTON  SINGLE-PARTICLE WF NOT FOUND ON THE RECORD FILE",1X,10(1H/))')
        Write(NFIPRI,'(1X,10(1H/)," RERUN THE BLOCKING CALCULATION WITH INIBLP=1           ",1X,10(1H/))')
        Stop ' PROTON SINGLE-PARTICLE WF NOT FOUND IN RECORD'
     End If

   End Subroutine read_hfodd

   !---------------------------------------------------------------------!
   !  Subroutine reads lists all the files in 'directory_name', extracts !
   !  all files of the form 'HFODD_XXXXXXXX.REC' (with XXXXXXXX a number)!
   !  and checks that all files with the number between 1 and N are      !
   !  present. It then broadcasts a boolean flag indicating if all files !
   !  are there, and a vector of integers giving the breakdown.          !
   !---------------------------------------------------------------------!
   Subroutine check_directory(directory_name,N)

      Integer, INTENT(IN) :: N
      Character(Len=132), INTENT(IN) :: directory_name

      Integer, Parameter :: ndatin=85
      Integer :: i,exitstat,iocheck,ierr,mpi_rank,mpi_err,numero,count
      Character(Len=80) :: filein
      Character(Len=150) :: commande
      Character(Len=132) :: string

      mpi_rank=0; mpi_err=0

      ! Getting the rank of the current process
      Call mpi_comm_rank(MPI_COMM_WORLD, mpi_rank, mpi_err)

      If(mpi_rank.Eq.0) Then

     	 commande = 'ls '//Trim(directory_name)// ' > toto'; exitstat=0
      	 Call system(commande)
      	 If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," - Executed command: ",a)') &
      	                           mpi_rank,commande

      	 ! open file
      	 Open(ndatin,file='toto',status='old',form='formatted',iostat=ierr)
         If(ierr.Ne.0) Then
            Write(6,'("Error in opening file : toto")')
            Stop 'Error in check_directory - I/O'
         End if

         ! Read file, extract number of file if relevant, and fill table
         iocheck=0
         Do While(iocheck.Eq.0)
            Read(ndatin,*,Iostat=iocheck) string
            numero = 0
            Call get_numero(string,numero)
            If(debug.Ge.2) Write(6,'("numero = ",i10)') numero
            If(numero.Gt.0) table_fichiers(numero)=1
         End Do

      	 Close(ndatin)
      	 If(debug.Ge.1) Write(6,'("Done reading with the file")')

         filled=.True.; count=0
         Do i=1,N
            If(table_fichiers(i).Eq.0) Then
               count=count+1
            End If
         End Do
         If(count.Gt.0) filled=.False.

      	 If(debug.Ge.2) Then
      	    If(filled) Then
               Write(6,'("All files are there")')
            Else
               Write(6,'(i10," files are missing")') count
            End If
         End If

      End If

      ! Broadcast flag telling if array is filled and table of missing files
      Call MPI_Bcast(filled,1,MPI_LOGICAL,0,MPI_COMM_WORLD,mpi_err)
      Call MPI_Bcast(table_fichiers,N,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err)

      ! Wait for all threads
      Call mpi_barrier(MPI_COMM_WORLD,mpi_err)

    End Subroutine check_directory

   !---------------------------------------------------------------------!
   !  Subroutine returns in 'numero' the file number from a string of    !
   !   the form 'HFODD_XXXXXXXX.REC' (with XXXXXXXX the number).         !
   !---------------------------------------------------------------------!
    Subroutine get_numero(string,numero)

      Character(Len=132), INTENT(IN) :: string
      Integer, INTENT(INOUT) :: numero
      Integer :: start
      Character(Len=4) :: dum_4
      Character(Len=6) :: dum_6

      start = Index(string, 'HFODD_', .False.)
      If (start.Ge.0) Then
          Read(string(start:start+18), '(A6,I8.8,A4)') dum_6,numero,dum_4
      End If

    End Subroutine get_numero

   !---------------------------------------------------------------------!
   !  Subroutine copies all files that have not been copied yet.         !
   !---------------------------------------------------------------------!
    Subroutine fill_holes(N,IF_THO,LIPKIN,LIPKIP)

      Integer, INTENT(IN) :: N,IF_THO,LIPKIN,LIPKIP

      Integer :: jj,exitstat,mpi_rank,mpi_world,mpi_err
      Character(Len=150) :: commande

      mpi_rank=0; mpi_world=1; mpi_err=0

      ! Getting the rank of the current process
      Call mpi_comm_rank(MPI_COMM_WORLD, mpi_rank, mpi_err)
      Call mpi_comm_size(MPI_COMM_WORLD, mpi_world,mpi_err)

      Do jj=1,N

      	 If(Mod(mpi_rank,mpi_world).Eq.Mod(jj-1,mpi_world)) Then

            If(table_fichiers(jj).Eq.0) Then

               commande = 'cp '//fichier_rec_old(jj)//' '//fichier_rec_new(jj); exitstat=0
               Call system(commande)
               If(debug.Ge.2) Write(6,'("fill_holes: mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
                                                     mpi_rank,jj,commande
               If(IF_THO.Ge.1) Then
                  commande = 'cp '//fichier_tho_old(jj)//' '//fichier_tho_new(jj); exitstat=0
                  Call system(commande)
                  If(debug.Ge.2) Write(6,'("fill_holes: mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
                                                        mpi_rank,jj,commande
               End If
               If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then
                  commande = 'cp '//fichier_lic_old(jj)//' '//fichier_lic_new(jj); exitstat=0
                  Call system(commande)
                  If(debug.Ge.2) Write(6,'("fill_holes: mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') &
                                                        mpi_rank,jj,commande
               End If

            End If

         End If

      End Do

      ! Wait for all threads
      Call mpi_barrier(MPI_COMM_WORLD,mpi_err)

   End Subroutine fill_holes

 End Module hfodd_mpimanager
