xref: /phasta/M2N/src/readnblk.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)

!  readnblk.f (pronounce "Reed and Block Dot Eff") contains:
!
!    module readarrays ("Red Arrays") -- contains the arrays that
!     are read in from binary files but not immediately blocked
!     through pointers.
!
!    subroutine readnblk ("Reed and Block") -- allocates space for
!     and reads data to be contained in module readarrays.  Reads
!     all remaining data and blocks them with pointers.
!


      module readarrays

      real*8, allocatable :: point2x(:,:)
      real*8, allocatable :: qold(:,:)
      real*8, allocatable :: dwal(:)
      real*8, allocatable :: errors(:,:)
      real*8, allocatable :: ybar(:,:)
      real*8, allocatable :: yphbar(:,:,:)
      real*8, allocatable :: vort(:,:)
      real*8, allocatable :: uold(:,:)
      real*8, allocatable :: acold(:,:)
      integer, allocatable :: iBCtmp(:)
      real*8, allocatable :: BCinp(:,:)

      integer, allocatable :: point2ilwork(:)
      integer, allocatable :: nBC(:)
      integer, allocatable :: point2iper(:)
      integer, allocatable :: point2ifath(:)
      integer, allocatable :: point2nsons(:)

      end module



      subroutine readnblk
!
      use readarrays
      include "commonM2N.h"
      include "mpif.h"

!
      real*8, allocatable :: xread(:,:), qread(:,:), acread(:,:)
      real*8, allocatable ::  qreadN(:), qoldN(:,:)
      real*8, allocatable :: pID(:), vID(:), dwalN(:), ybarN(:,:)
      real*8, allocatable :: errorsN(:,:), yphbarN(:,:,:), vortN(:,:)
      real*8, allocatable :: qTota2as(:,:), qTota2ar(:,:)
      real*8, allocatable :: qrecv(:), qsend(:)
      real*8 :: qreadN1
      integer, allocatable :: iperread(:), iBCtmpread(:)
      integer, allocatable :: ilworkread(:), nBCread(:)
      integer, allocatable :: indexpart(:)
      integer, allocatable :: getinfo(:)
      character*10 cname2
      character*8 mach2
      character*30 fmt1
      character*255 fname1,fnamer, fnamelr
      character*255 warning
      integer :: igeomBAK, ibndc, irstin, ierr
      integer ::  icountN(numpe) ! integers read from headers
      integer :: intfromfile(50) ! integers read from headers
      integer :: maxvID(numpe),maxvIDout(numpe)
      integer :: nmap, ndim
      integer :: maxicountN, maxicountNglob, numvar, sumvar
      integer :: ncount, ifill, tag
      integer :: stat (MPI_STATUS_SIZE)
      integer :: ndofdwal,ndofybar,ndoferrors,ndofyphbar,ndofvort
      logical :: exmap, exinput
      integer ::  my_local_comm, my_local_size, my_local_rank
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cccccccccccccccccccccccccccccccccccccccc New PhastaIO Definition Part ccccccccccccccccccccccccccccccccccccccccc

      integer :: descriptor, descriptorG, GPID, color, nfiles, nfields
      integer :: numparts, nppf, islesseqN
      integer :: ierr_io, numprocs, itmp, itmp2
      integer :: i, nonzero
      integer :: ireducemethod, isbinary, irstart, irstartmap, iybar
      integer :: ierror, numphavg, ivort, idwal, idebug, iphavg, irankN
      character*255 fname2, temp2
      character*64 temp1

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

!
!  Set some parameters for the reduce operation
!
!     0 for communication based approach with mpi_send and mpi_recv (recommended), 1 for file based approach (not maintained any more)
      ireducemethod = 0

!     1 to write the restartMap* files in binary format, 0 in ascii (in conjunction with ireducemethod = 1)
      isbinary = 1
!
!     Read the input parameter for the the M2N reduction
!
!      open(unit=72,file='numstart.dat',status='old')
!      read(72,*) irstart
!      close(72)

      if(myrank == 0) then
        write(*,*) 'Reading M2N_input.dat'
        fnamer='M2N_input.dat'
        inquire(file=fnamer,exist=exinput)
        if(exinput) then
          open(unit=72,file=fnamer,status='old')
          read(72,*) irstart
          read(72,*) irstartmap
          read(72,*) iybar
          read(72,*) ierror
          read(72,*) numphavg
          read(72,*) ivort
          read(72,*) idwal
          read(72,*) idebug
          close(72)
        else
           write(*,*) 'ERROR: Input file ',
     &                 trim(fnamer),' does not exist!'
        endif

      endif

      call mpi_barrier(mpi_comm_world, ierr)
      call mpi_bcast(exinput,1,MPI_LOGICAL,0,mpi_comm_world,ierr)
      if(.not. exinput) then ! M2N_input.dat does not exist. Quit
        return
      else ! broadcast the information read by rank 0
        call mpi_bcast(irstart,1,MPI_INTEGER,0,mpi_comm_world,ierr)
        call mpi_bcast(irstartmap,1,MPI_INTEGER,0,mpi_comm_world,ierr)
        call mpi_bcast(iybar,1,MPI_INTEGER,0,mpi_comm_world,ierr)
        call mpi_bcast(ierror,1,MPI_INTEGER,0,mpi_comm_world,ierr)
        call mpi_bcast(numphavg,1,MPI_INTEGER,0,mpi_comm_world,ierr)
        call mpi_bcast(ivort,1,MPI_INTEGER,0,mpi_comm_world,ierr)
        call mpi_bcast(idwal,1,MPI_INTEGER,0,mpi_comm_world,ierr)
        call mpi_bcast(idebug,1,MPI_INTEGER,0,mpi_comm_world,ierr)
      endif

      if(myrank == 0 ) then
      ! Print some info
        write(*,*) 'The solution field is reduced by default'

        if(iybar .gt. 0) then
          write(*,*) 'The ybar field will also be reduced'
          iybar = 1 ! security
        else
          write(*,*) 'The ybar field will NOT be reduced'
        endif

        if(ierror .gt. 0) then
          write(*,*) 'The error field will also be reduced'
          ierror = 1 ! security
        else
          write(*,*) 'The error field will NOT be reduced'
        endif

        if(numphavg .gt. 0) then
          write(*,*) 'The phase average fields (',numphavg,
     &                ') will also be reduced'
        else
          write(*,*) 'The phase average fields will NOT be reduced'
        endif

        if(ivort .gt. 0) then
          write(*,*) 'The vorticity field will also be reduced'
          ivort = 1 ! security
        else
          write(*,*) 'The vorticity field will NOT be reduced'
        endif

        if(idwal .gt. 0) then
          write(*,*) 'The dwal field will also be reduced'
          idwal = 1 ! security
        else
          write(*,*) 'The dwal field will NOT be reduced'
        endif
        write(*,*) ''
      endif

      call mpi_barrier(mpi_comm_world, ierr)
      if(myrank == 0) then
        write(*,*) 'Reading the geombc-dat files'
      endif

      lstep=irstart ! in case restart files have no fields

      nfiles = nsynciofiles
      numparts = numpe !This is the common settings. Beware if you try to compute several parts per process

      color = int(myrank/(numparts/nfiles)) !Should call the color routine in SyncIO here
      itmp2 = int(log10(float(color+1)))+1
      write (temp2,"('(''geombc-dat.'',i',i1,')')") itmp2
      write (fnamer,temp2) (color+1)
      fnamer = trim(fnamer)//char(0)

      itwo=2
      ione=1
      ieleven=11
      itmp = int(log10(float(myrank+1)))+1

      call queryphmpiio(fnamer, nfields, nppf);
      if (myrank == 0) then
        write(*,*) 'Number of fields in geombc-dat: ',nfields
        write(*,*) 'Number of parts per file geombc-dat: ',nppf
      endif
      call initphmpiio( nfields, nppf, nfiles, igeom,
     & 'read'//char(0))
      call openfile( fnamer, 'read'//char(0), igeom )

      write (temp1,"('(''number of nodes@'',i',i1,',A1)')") itmp
      write (fname2,temp1) (myrank+1),'?'
      call readheader(igeom,fname2//char(0),numnp,ione,
     & 'integer'//char(0), iotype)

      write (temp1,"('(''number of modes@'',i',i1,',A1)')") itmp
      write (fname2,temp1) (myrank+1),'?'
      call readheader(igeom,fname2//char(0),nshg,ione,
     & 'integer'//char(0), iotype)

      write (temp1,"('(''number of interior elements@'',i',i1,',A1)')")
     &       itmp
      write (fname2,temp1) (myrank+1),'?'
      call readheader(igeom,fname2//char(0),numel,ione,
     & 'integer'//char(0), iotype)

      write (temp1,"('(''number of boundary elements@'',i',i1,',A1)')")
     &       itmp
      write (fname2,temp1) (myrank+1),'?'
      call readheader(igeom,fname2//char(0),numelb,ione,
     & 'integer'//char(0),iotype)

      write (temp1,
     & "('(''maximum number of element nodes@'',i',i1,',A1)')") itmp
      write (fname2,temp1) (myrank+1),'?'
      call readheader(igeom,fname2//char(0),nen,ione,
     &'integer'//char(0),iotype)

      write (temp1,"('(''number of interior tpblocks@'',i',i1,',A1)')")
     &       itmp
      write (fname2,temp1) (myrank+1),'?'
      call readheader(igeom,fname2//char(0),nelblk,ione,
     & 'integer'//char(0) ,iotype)

      write (temp1,"('(''number of boundary tpblocks@'',i',i1,',A1)')")
     &       itmp
      write (fname2,temp1) (myrank+1),'?'
      call readheader(igeom,fname2//char(0),nelblb,ione,
     & 'integer'//char(0), iotype)

      write (temp1,
     & "('(''number of nodes with Dirichlet BCs@'',i',i1,',A1)')") itmp
      write (fname2,temp1) (myrank+1),'?'
      call readheader(igeom,fname2//char(0),numpbc,ione,
     & 'integer'//char(0),iotype)

      write (temp1,"('(''number of shape functions@'',i',i1,',A1)')")
     &       itmp
      write (fname2,temp1) (myrank+1),'?'
      call readheader(igeom,fname2//char(0),ntopsh,ione,
     & 'integer'//char(0),iotype)

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
!.... calculate the maximum number of boundary element nodes
!
      nenb = 0
      do i = 1, melCat
         if (nen .eq. nenCat(i,nsd)) nenb = max(nenCat(i,nsd-1), nenb)
      enddo
!
!
!.... setup some useful constants
!
      I3nsd  = nsd / 3          ! nsd=3 integer flag
      E3nsd  = float(I3nsd)     ! nsd=3 real    flag
!
      if(matflg(1,1).lt.0) then
         nflow = nsd + 1
      else
         nflow = nsd + 2
      endif
      ndof   = nsd + 2
      nsclr=impl(1)/100
      ndof=ndof+nsclr           ! number of sclr transport equations to solve

      ndofBC = ndof + I3nsd     ! dimension of BC array
      ndiBCB = 2                ! dimension of iBCB array
      ndBCB  = ndof + 1         ! dimension of BCB array
!
      nsymdf = (ndof*(ndof + 1)) / 2 ! symm. d.o.f.'s
!
!.... ----------------------> Communication tasks <--------------------
!
      call closefile( igeom, "read"//char(0) )
      call finalizephmpiio( igeom )

!
!.... Read restart files for the solution, dwal, error and ybar
!
      call mpi_barrier(mpi_comm_world, ierr)
      if(myrank == 0) then
        write(*,*)'Reading the restart-dat files for the s-d-e-y fields'
      endif

      itmp=1
      if (irstart .gt. 0) itmp = int(log10(float(irstart+1)))+1

      write (fmt1,"('(''restart-dat.'',i',i1,',1x)')") itmp

      write (fnamer,fmt1) irstart
      fnamer = trim(fnamer) // cname2(color+1)

      call queryphmpiio(fnamer//char(0), nfields, nppf);
      if (myrank == 0) then
        write(*,*) 'Number of fields in restart-dat: ',nfields
        write(*,*) 'Number of parts per file restart-dat: ',nppf
      endif
      call initphmpiio(nfields,nppf,nfiles,descriptor,
     & 'read'//char(0))
      call openfile( fnamer//char(0) ,
     & 'read'//char(0), descriptor )

      ithree=3
!      call creadlist(irstin,ithree,nshg2,ndof2,lstep)

      itmp = int(log10(float(myrank+1)))+1
      write (temp1,"('(''solution@'',i',i1,',A1)')") itmp
      write (fname1,temp1) (myrank+1),'?'
      fname1 = trim(fname1)

!      print *, "Solution is : ", fname1

      intfromfile=0
      call readheader(descriptor,fname1//char(0) ,intfromfile,
     & ithree,'integer'//char(0), iotype)
!
!.... read the values of primitive variables into q
!

!      print *, "intfromfile(1) is ", intfromfile(1)
!      print *, "intfromfile(2) is ", intfromfile(2)
!      print *, "intfromfile(3) is ", intfromfile(3)

      if(intfromfile(1).ne.0) then
         nshg2=intfromfile(1)
         ndof2=intfromfile(2)
         lstep=intfromfile(3)
         if(nshg2.ne.nshg) then
           write(*,*) 'ERROR: mistmatch between nshg '//
     &      'from the geombc and restart files',nshg,nshg2
         endif

         ndof=ndof2
         allocate( qold(nshg,ndof2) )
         allocate( qread(nshg,ndof2) )
         iqsiz=nshg*ndof2
         call readdatablock(descriptor,fname1//char(0),qread,iqsiz,
     &                         'double'//char(0),iotype)
         qold(:,1:ndof)=qread(:,1:ndof)
         deallocate(qread)
      else
         if (myrank.eq.master) then
            if (matflg(1,1).eq.0) then ! compressible
               warning='WARNING: Solution is set to zero '//
     &                   '(with p and T to one)'
            else
               warning='WARNING: Solution is set to zero'
            endif
            write(*,*) warning
         endif
         qold=zero
         if (matflg(1,1).eq.0) then ! compressible
            qold(:,1)=one ! avoid zero pressure
            qold(:,nflow)=one ! avoid zero temperature
         endif
      endif

!
!  Above this line is the usual loading of fields in readnblk.  Quite a bit of stuff has been removed.
!


!
! follow the usual convention for loading the ybar
!
      ndofybar=0
      if(iybar == 1) then
        itmp = int(log10(float(myrank+1)))+1
        write (temp1,"('(''ybar@'',i',i1,',a1)')") itmp
        write (fname1,temp1) (myrank+1),'?'
        fname1 = trim(fname1)
        intfromfile=0
        call readheader(descriptor,fname1//char(0),intfromfile,
     &     ithree,'integer'//char(0),iotype)

        !nshg2=intfromfile(1)
        ndofybar=intfromfile(2)
        !lstep=intfromfile(3)

        if(ndofybar .ne. 0) then
          allocate( ybar(nshg,ndofybar) )
          iqsiz=nshg*ndofybar
          call readdatablock(descriptor,fname1//char(0) ,ybar,iqsiz,
     &                   'double'//char(0),iotype)
        else
          if(myrank==0) then
            write(*,*) 'WARNING: ybar is missing in the restart files'
          endif
          iybar = 0
        endif
      endif

!
! follow the usual convention for loading the error
!
      ndoferrors=0
      if(ierror == 1) then
        itmp = int(log10(float(myrank+1)))+1
        write (temp1,"('(''errors@'',i',i1,',a1)')") itmp
        write (fname1,temp1) (myrank+1),'?'
        fname1 = trim(fname1)
        intfromfile=0
        call readheader(descriptor,fname1//char(0),intfromfile,
     &     ithree,'integer'//char(0),iotype)

        !nshg2=intfromfile(1)
        ndoferrors=intfromfile(2)
        !lstep=intfromfile(3)

        if(ndoferrors .ne. 0) then
          allocate( errors(nshg,ndoferrors) )
          iqsiz=nshg*ndoferrors
          call readdatablock(descriptor,fname1//char(0),errors,iqsiz,
     &                   'double'//char(0),iotype)
        else
          if(myrank==0) then
            write(*,*) 'WARNING: errors is missing in the restart files'
          endif
          ierror = 0
        endif
      endif

!
!     Read the phase_average fields
!
      ndofyphbar=0
      if(numphavg .gt. 0) then
        do iphavg = 1,numphavg
          itmp = int(log10(float(myrank+1)))+1
          itmp2 = int(log10(float(iphavg)))+1
          write (temp1,
     &            "('(''phase_average'',i',i1,',''@'',i',i1,',A1)')")
     &             itmp2, itmp
          write (fname1,temp1) iphavg,(myrank+1),'?'
          fname1 = trim(fname1)
          intfromfile=0
          call readheader(descriptor,fname1//char(0),intfromfile,
     &                  ithree,'integer'//char(0),iotype)

          !nshg=intfromfile(1) !Do not use nshg and ndof from common.h here!
          ndofyphbar=intfromfile(2)
          !lstep=intfromfile(3)

          if(ndofyphbar.ne.0) then
            ! Allocate some memory for the first ts only
            if(iphavg==1) then
              allocate( yphbar(nshg,ndofyphbar,numphavg) )
            endif

            allocate( qread(nshg,ndofyphbar)  ) ; qread = 0.d0
            iqsiz = nshg*ndofyphbar
            call readdatablock(descriptor,fname1//char(0),qread,iqsiz,
     &                     'double'//char(0),iotype)
            yphbar(:,:,iphavg) = qread(:,:)
            deallocate(qread)
          else
            if(myrank==0) then
              write(*,*)'WARNING: phase_average is missing '//
     &            'in the restart files'
            endif
            numphavg = 0
            if(iphavg > 1) then
              deallocate(yphbar)
            endif
            exit
          endif
        enddo
      endif

!
! follow the usual convention for loading the vorticity field
!
      ndofvort=0
      if(ivort == 1) then
        itmp = int(log10(float(myrank+1)))+1
        write (temp1,"('(''vorticity@'',i',i1,',a1)')") itmp
        write (fname1,temp1) (myrank+1),'?'
        fname1 = trim(fname1)
        intfromfile=0
        call readheader(descriptor,fname1//char(0),intfromfile,
     &     ithree,'integer'//char(0),iotype)

        !nshg2=intfromfile(1)
        ndofvort=intfromfile(2)
        !lstep=intfromfile(3)

        if(ndofvort .ne. 0) then
          allocate( vort(nshg,ndofvort) )
          iqsiz=nshg*ndofvort
          call readdatablock(descriptor,fname1//char(0) ,vort ,iqsiz,
     &                   'double'//char(0),iotype)
        else
          if(myrank==0) then
            write(*,*) 'WARNING: vorticity is missing '//
     &                 'in the restart files'
          endif
          ivort = 0
        endif
      endif


!
! follow the usual convention for loading the d2wal
!
      ndofdwal=0
      if(idwal == 1 ) then
        write (temp1,"('(''dwal@'',i',i1,',a1)')")
     &       itmp
        write (fname1,temp1) (myrank+1),'?'
        fname1 = trim(fname1)
        intfromfile=0
        call readheader(descriptor,fname1//char(0),intfromfile,
     &     ithree,'integer'//char(0),iotype)

         !nshg2=intfromfile(1)
         ndofdwal=intfromfile(2)
         !lstep=intfromfile(3)

        if(ndofdwal .ne. 0) then
          if(ndofdwal.ne.1) then
            warning='WARNING: ndofdwal not equal 1'
            write(*,*) warning, ndofdwal
          endif
          allocate( dwal(nshg) )
          iqsiz=nshg
          call readdatablock(descriptor,fname1//char(0), dwal, iqsiz,
     &                   'double'//char(0),iotype)
        else
          if(myrank==0) then
            write(*,*) 'WARNING: dwal is missing in the restart files'
          endif
          idwal = 0
        endif
      endif

!
!.... close c-binary files
!
      call closefile( descriptor, "read"//char(0) )
      call finalizephmpiio( descriptor )

!
!.... Read restart files for the vtx mapping between M and N parts
!
      call mpi_barrier(mpi_comm_world, ierr)
      if(myrank == 0) then
        write(*,*) 'Reading the restart-dat files for vtx mapping'
      endif

      itmp=1
      if (irstartmap .gt. 0) itmp = int(log10(float(irstartmap+1)))+1

      write (fmt1,"('(''restart-dat.'',i',i1,',1x)')") itmp

      write (fnamer,fmt1) irstartmap
      fnamer = trim(fnamer) // cname2(color+1)

      call queryphmpiio(fnamer//char(0), nfields, nppf);
      if (myrank == 0) then
        write(*,*) 'Number of fields in restart-dat: ',nfields
        write(*,*) 'Number of parts per file restart-dat: ',nppf
      endif
      call initphmpiio(nfields,nppf,nfiles,descriptor,
     & 'read'//char(0))
      call openfile( fnamer//char(0) ,
     & 'read'//char(0), descriptor )

!
! follow the usual convention for loading the partid
!
      itmp = int(log10(float(myrank+1)))+1
       write (temp1,"('(''mapping_partid@'',i',i1,',A1)')")
     &          itmp
       write (fname1,temp1) (myrank+1),'?'
       fname1 = trim(fname1)
       intfromfile=0
       call readheader(descriptor,fname1//char(0),intfromfile,
     &        ithree,'integer'//char(0),iotype)

       nshg2=intfromfile(1)
       nmap=intfromfile(2)
       !lstep=intfromfile(3) !lstep comes from another restart file
       if(nmap.ne.1) then
          warning='WARNING mapping_partid not equal 1'
          write(*,*) warning , nmap, intfromfile(1), intfromfile(2)
       endif

       allocate( pID(nshg) )
       iqsiz=nshg
       call readdatablock(descriptor,fname1//char(0) ,pID,iqsiz,
     &                   'double'//char(0),iotype)

!
! follow the usual convention for loading the vtxid
!
        write (temp1,"('(''mapping_vtxid@'',i',i1,',a1)')")
     &           itmp
        write (fname1,temp1) (myrank+1),'?'
        fname1 = trim(fname1)
        intfromfile=0
        call readheader(descriptor,fname1//char(0),intfromfile,
     &     ithree,'integer'//char(0),iotype)

       nshg2=intfromfile(1)
       nmap=intfromfile(2)
       !lstep=intfromfile(3) !lstep comes from another restart file
       if(nmap.ne.1) then
          warning='warning mapping_partid not equal 1'
          write(*,*) warning , nmap
       endif

       allocate( vid(nshg) )
       iqsiz=nshg
       call readdatablock(descriptor,fname1//char(0) ,vid,iqsiz,
     &                   'double'//char(0),iotype)

!
!.... close c-binary files
!
      call closefile( descriptor, "read"//char(0) )
      call finalizephmpiio( descriptor )

!
!  Debugging information
!
!      minvid=minval(vid)
!      write(*,*) minvid,myrank
!      do i=1,nshg
!      write(800+myrank,*) pID(i),vID(i),qold(i,1), qold(i,2)
!      ivID=vID(i)+0.1
!      ipID=pID(i)+0.1
!      write(900+myrank,*) ipID,ivID,qold(i,1), qold(i,2)
!      enddo

!
! here starts the mapping code
!

!  first, on each of the M parts, count the number of vertices on each of the
!  N parts and put the result in icountN for each of the M parts
!
      call mpi_barrier(mpi_comm_world, ierr)
      if(myrank == 0) then
        write(*,*) 'Computing the size of the mapped partition'
      endif

      maxPID=0
      icountN=0
      do i=1,nshg
          ifileN=pID(i)+1.1
          icountN(ifileN)=icountN(ifileN)+1
          maxPID=max(maxPID,ifileN)
      enddo
!
!     ALLREDUCE so that every processor knows N the size of the mapped partition
!
      call MPI_ALLREDUCE (maxPID, maxPIDout, 1,
     &          MPI_INT, MPI_MAX, mpi_comm_world, ierr)

      irankN=maxPIDout

      call mpi_barrier(mpi_comm_world, ierr)
      if(myrank.eq.0) then
        write (*,*) 'M2N reduction with M =',
     &                        numpe, 'and N =',irankN
      endif
!
!     DEBUG echo to the screen the number of number of vtx that will be mapped to
!     each of the N for each of the M parts
!
      if (idebug == 1) then
        call mpi_barrier(mpi_comm_world, ierr)
        do i=1,numpe
          if(myrank == i) then
            write(*,926) myrank, (icountN(j),j=1,irankN)
          endif
          call mpi_barrier(mpi_comm_world, ierr)
        enddo
      endif

!
!     Compute the max vertex ID for every N part
!
      call mpi_barrier(mpi_comm_world, ierr)
      if(myrank == 0) then
        write(*,*) 'Computing the max vertex ID for every N part'
      endif

      maxvID(:)=0 !only maxvID(1:irankN) should have non-zero entries at the end of the program
      do i=1,nshg
        ifileN=pID(i)+1.1    !shift pID to start from 1
        intvID=(vID(i)+1.1)  !shift vID to start from 1
        maxvID(ifileN)=max(maxvID(ifileN),intvID)
      enddo

!
!     Sanity check for maxvID
!
      do i=irankN+1,numpe
        if (maxvID(i) .ne. 0) then
           write(*,*) 'ERROR: maxvID(',i,') not equal to 0 on rank',
     &                 myrank
        endif
      enddo
!
!     ALLREDUCE so that every processor knows the number of vtx on each of the N
!     parts in the reduced partition. No need to reduce maxvID(irankN+1:numpe), as it should be 0
!
      call MPI_ALLREDUCE (maxvID, maxvIDout, irankN,
     &          MPI_INT, MPI_MAX, mpi_comm_world, ierr)

!
!    echo debug information to the screen
!
      if (idebug == 1) then
        if(myrank.eq.0) write(*,*) 'Writing maxvID on each of N ranks'
        if(myrank.eq.0) write(*,926) (maxvIDout(j),j=1,irankN)
      endif

!     Beware! maxvIDout does not represent the largest vID present in the parts
!     but the largest vID referenced in the parts through the mapping.
!     On part boundaries, it is possible to have unreferenced vID that will need
!     to be fixed through P2P communication.

!
!     Allocate the memory for qoldN on N parts handled by ranks 1..N.
!     Memory available per comupte node could be improved by chosing others ranks
!
      call mpi_barrier(mpi_comm_world, ierr)
      if(myrank == 0) then
        write(*,*) 'Allocating memory for the reduced s-d-e-y vectors'
      endif

      if(myrank.lt.irankN) then
        nshgN=maxvIDout(myrank+1) ! number of vtx on my part computed above

        allocate(qoldN(nshgN,ndof))! size of my final mapped field array
        qoldN=-9.87654321e32       ! set to absurd, easily recognized value to find holes

        if(iybar == 1) then
          allocate(ybarN(nshgN,ndofybar))   ! size of my final mapped field array
          ybarN=-9.87654321e32       ! set to absurd, easily recognized value to find holes
        endif

        if(ierror == 1) then
          allocate(errorsN(nshgN,ndoferrors))   ! size of my final mapped field array
          errorsN=-9.87654321e32       ! set to absurd, easily recognized value to find holes
        endif

        if(numphavg .gt. 0) then
          allocate(yphbarN(nshgN,ndofyphbar,numphavg))   ! size of my final mapped field array
          yphbarN=-9.87654321e32       ! set to absurd, easily recognized value to find holes
        endif

        if(ivort == 1) then
          allocate( vortN(nshgN,ndofvort))
          vortN=-9.87654321e32
        endif

        if(idwal == 1) then
          allocate(dwalN(nshgN))   ! size of my final mapped field array
          dwalN=-9.87654321e32       ! set to absurd, easily recognized value to find holes
        endif


      else ! Not used but sane, as passed to write_M2N
        allocate(qoldN(1,1))
        if(iybar == 1) then
          allocate(ybarN(1,1))
        endif
        if(ierror == 1) then
          allocate(errorsN(1,1))
        endif
        if(numphavg .gt. 0) then
          allocate(yphbarN(1,1,1))
        endif
        if(ivort == 1) then
          allocate(vortN(1,1))
        endif
        if(idwal == 1) then
          allocate(dwalN(1))
        endif
      endif

!
!     Most important routine: build the qoldN array for the reduction of the solution from M to N
!
      call mpi_barrier(mpi_comm_world, ierr)
      if(myrank == 0) then
        write(*,*) 'Building the reduced field vectors'
      endif

      if (ireducemethod  == 0) then
!
!       Communication approach through mpi_isend and mpi_irecv to limit memory consumption
!       and avoid disk access
!
        numvar = 1+ndof+ndofybar+ndoferrors
     &           +(numphavg*ndofyphbar)+ndofvort+ndofdwal
!        write(*,*) numvar,ndof,ndofybar,ndoferrors,numphavg,ndofyphbar,
!     &             ndofvort
        call mpi_barrier(mpi_comm_world, ierr)
        !1 for vtxid, ndof for solution, ndofdwal for dwal, ndoferrors for error, ndofybar for ybar

        do irN = 0,irankN-1  ! every irN rank 1..N receives information from


          ! Get which of the numpe-1 ranks have some information to communicate to rank irN
          if (myrank == irN) then
            allocate(getinfo(numpe))
            getinfo(:) = 0
          endif

          ifill = icountN(irN+1)
          call mpi_gather(ifill, 1, MPI_INT, getinfo, 1, MPI_INT,
     &                   irN, mpi_comm_world, ierr)

!          if (myrank == irN) then
!            write(*,*) getinfo(:)
!          endif

          do irM = 0,numpe-1 ! all the irM ranks 1..M

! DEBUG
!            if (myrank == 0) then
!              write(*,*) 'NEW data exchange between irN:',
!     &                  irN,'and irM:',irM
!            endif
!            call mpi_barrier(mpi_comm_world, ierr)
!DEBUG

            if (myrank == irN .and. myrank == irM) then ! no need to send the data. Already there so get them
!             write(*,*) '  myrank=irN=irM=',irN,' so treat data locally'
              do i=1,nshg
                ifileN = pID(i)+1.1 !shift to start from 1
                if(ifileN == irN+1) then ! these data need to be packed because rank irN is expecting them
                  intvID = vID(i)+1.1 ! shift to start from 1
                  ! Solution qoldN
                  do idof=1,ndof
                    qoldN(intvID,idof) = qold(i,idof)
                  enddo
                  ! ybarN
                  if(iybar==1) then
                    do idof=1,ndofybar
                      ybarN(intvID,idof) = ybar(i,idof)
                    enddo
                  endif
                  ! errorsN
                  if(ierror==1) then
                    do idof=1,ndoferrors
                      errorsN(intvID,idof) = errors(i,idof)
                    enddo
                  endif
                  ! phase average
                  if(numphavg .gt. 0) then
                    do iphavg = 1,numphavg
                      do idof=1,ndofyphbar
                        yphbarN(intvID,idof,iphavg) =
     &                        yphbar(i,idof,iphavg)
                      enddo
                    enddo
                  endif
                  ! vortN
                  if(ivort==1) then
                    do idof=1,ndofvort
                      vortN(intvID,idof) = vort(i,idof)
                    enddo
                  endif
                 ! Distance to the wall dwalN
                  if(idwal==1) then
                    dwalN(intvID) = dwal(i)
                  endif
               endif
              enddo

            else !myrank ne to irN and irM
              tag=irN
              if(myrank == irN) then
                ! rank irN receives data from rank j and update qoldN.
!                call mpi_recv(ifill, 1, MPI_INTEGER, irM, tag,
!     &                                    mpi_comm_world, stat, ierr)
!                write(*,*) '  ', irN,'receives',ifill,getinfo(irM+1),
!     &                     'data from',irM
                ifill = getinfo(irM+1)
                if(ifill .gt. 0) then
                  ! there are really data to receive
!                  write(*,*) '  ', irN,'receives',ifill,'data from',irM
                  allocate(qrecv(ifill*numvar))
                  qrecv = -9.87654321e32
                  call mpi_recv(qrecv, ifill*numvar, MPI_DOUBLE, irM,
     &                                tag, mpi_comm_world, stat, ierr)
                  do i=1,ifill
                    intvID = qrecv( (i-1)*numvar+1 )+1.1 ! shift to start from 1
                    ! Solution qold
                    do idof=1,ndof
                      sumvar = 1+idof
                      qoldN(intvID,idof) = qrecv( (i-1)*numvar+sumvar )
                    enddo
                    ! ybarN
                    if(iybar==1) then
                      do idof=1,ndofybar
                        sumvar = 1+ndof+idof
                        ybarN(intvID,idof) = qrecv( (i-1)*numvar+sumvar)
                      enddo
                    endif
                    ! errorN
                    if(ierror==1) then
                      do idof=1,ndoferrors
                        sumvar = 1+ndof+ndofybar+idof
                        errorsN(intvID,idof) =
     &                            qrecv( (i-1)*numvar+sumvar)
                      enddo
                    endif
                    ! yphbarN
                    if(numphavg .gt. 0) then
                      do iphavg = 1,numphavg
                        do idof=1,ndofyphbar
                          sumvar = 1+ndof+ndofybar+ndoferrors+
     &                                    (iphavg-1)*ndofyphbar+idof
                          yphbarN(intvID,idof,iphavg) =
     &                                     qrecv( (i-1)*numvar+sumvar)
                        enddo
                      enddo
                    endif
                    ! vortN
                    if(ivort==1) then
                      do idof=1,ndofvort
                        sumvar = 1+ndof+ndofybar+ndoferrors+
     &                                    numphavg*ndofyphbar+idof
                        vortN(intvID,idof) =
     &                            qrecv( (i-1)*numvar+sumvar)
                      enddo
                    endif
                    ! Distance to the wall dwalN
                    if(idwal==1) then
                      sumvar = 1+ndof+ndofybar+ndoferrors+
     &                         numphavg*ndofyphbar+ndofvort+1
                      dwalN(intvID) = qrecv ( (i-1)*numvar+sumvar)
                    endif
                  enddo
                  deallocate(qrecv)
                endif

              elseif(myrank == irM) then ! rank irM sends data to rank irN
                ! First build the data to send to rank i for qoldN
                ifill = 0
                if (icountN(irN+1) .gt. 0) then !there are data to transfer
                  allocate(qsend(icountN(irN+1)*numvar))
                  qsend(:) = -9.87654321e32
                  do i=1,nshg
                    ifileN = pID(i)+1.1 !shift to start from 1
                    if(ifileN == irN+1) then
                      ! these data need to be packed because rank irN is expecting them
                      ifill = ifill+1
                      ! vtx ID
                      qsend( (ifill-1)*numvar+1) = vID(i) ! no shift to start from 1 yet
                      ! Solution qold
                      do idof=1,ndof
                        sumvar = 1+idof
                        qsend( (ifill-1)*numvar+sumvar) = qold(i,idof)
                      enddo
                      ! ybar
                      if(iybar==1) then
                        do idof=1,ndofybar
                          sumvar = 1+ndof+idof
                          qsend( (ifill-1)*numvar+sumvar) = ybar(i,idof)
                        enddo
                      endif
                     ! errors
                      if(ierror==1) then
                        do idof=1,ndoferrors
                          sumvar = 1+ndof+ndofybar+idof
                          qsend( (ifill-1)*numvar+sumvar) =
     &                                                 errors(i,idof)
                        enddo
                      endif
                      ! yphbarN
                      if(numphavg .gt. 0) then
                        do iphavg = 1,numphavg
                          do idof=1,ndofyphbar
                            sumvar = 1+ndof+ndofybar+ndoferrors+
     &                                    (iphavg-1)*ndofyphbar+idof
                            qsend( (ifill-1)*numvar+sumvar) =
     &                                            yphbar(i,idof,iphavg)
                          enddo
                        enddo
                      endif
                      ! vortN
                      if(ivort==1) then
                        do idof=1,ndofvort
                          sumvar = 1+ndof+ndofybar+ndoferrors+
     &                                    numphavg*ndofyphbar+idof
                          qsend( (ifill-1)*numvar+sumvar) =
     &                                          vort(i,idof)
                        enddo
                      endif
                      ! Distance to the wall dwal
                      if(idwal==1) then
                        sumvar = 1+ndof+ndofybar+ndoferrors+
     &                           numphavg*ndofyphbar+ndofvort+1
                        qsend( (ifill-1)*numvar+sumvar) = dwal(i)
                      endif
                    endif
                  enddo
                endif
                if(ifill .ne. icountN(irN+1)) then
                  write(*,*) 'ERROR with data set:', irM, irN,
     &                                   ifill, icountN(irN+1)
                endif
!                write(*,*) '  ', irM,'sends', ifill,'data to',irN
!                call mpi_send(ifill, 1 ,MPI_INTEGER, irN, tag,
!     &                                          mpi_comm_world, ierr) ! Send the size of the data
                if(ifill .gt. 0) then !if the size of the data is >0, send it for good
!                  write(*,*) '  ', irM,'sends', ifill,'data to',irN
                  call mpi_send(qsend(1), ifill*numvar, MPI_DOUBLE, irN,
     &                                      tag, mpi_comm_world, ierr)
                  deallocate(qsend)
                endif

              endif !if(myrank==irN)
            endif !if (myrank == irN = irM)

          enddo !loop over irM

          call mpi_barrier(mpi_comm_world,ierr)
          if (myrank == irN) then
            nonzero = 0
            do i = 1,numpe
              if (getinfo(i) .gt. 0) then
                nonzero = nonzero +1
              endif
            enddo
            write(*,*)"Part ",irN,"out of N = ",irankN," updated from",
     &                nonzero, "parts"
            deallocate(getinfo)
          endif
        enddo !loop over irN

!
!       Do it without files through mpialltoall
!
!      maxicountN=maxval(icountN(:))
!      call MPI_ALLREDUCE (maxicountN, maxicountNglob, 1,
!     &          MPI_INT, MPI_MAX, mpi_comm_world, ierr)
!
!      numvar = 1+ndof ! vID + ndof from solution
!      write(*,*) 'maxicountN:',myrank, maxicountN, maxicountNglob, nshg
!      call mpi_barrier(mpi_comm_world, ierr)
!      allocate(qTota2as(numvar*maxicountNglob,numpe))
!      qTota2as(:,:) = -9.87654321e32

!      allocate(indexpart(numpe)); indexpart(:)=0
!      do i=1,nshg
!        ifileN = pID(i)+1.1                              !shift pID to start from 1
!        indexpart(ifileN) = indexpart(ifileN)+1
!        qTota2as( (indexpart(ifileN)-1)*numvar+1 , ifileN ) = vID(i)   ! vID not yet shifted to start from 1 here
!        do j=1,ndof
!
!         if((indexpart(ifileN)-1)*numvar+1+j>numvar*maxicountNglob) then
!            write(*,*) 'ERROR DIMENSION1:',myrank,
!     &       ifileN, indexpart(ifileN),
!     &       (indexpart(ifileN)-1)*numvar+1+j, numvar*maxicountNglob
!          endif
!          if (ifileN > numpe ) then
!              write(*,*) 'ERROR DIMENSION2:',myrank, ifileN, numpe
!          endif
!          qTota2as( (indexpart(ifileN)-1)*numvar+1+j, ifileN )=qold(i,j)
!        enddo
!      enddo

!      allocate(qTota2ar(numvar*maxicountNglob,numpe))
!      qTota2ar(:,:) = -9.87654321e32

!      ncount = numvar*maxicountNglob
!      call mpi_alltoall(qTota2as(1,1), ncount, MPI_DOUBLE,
!     &                  qTota2ar(1,1), ncount, MPI_DOUBLE,
!     &                  mpi_comm_world, ierr)
!      if (ierr .ne. 0) then
!        write(*,*) 'Error with mpi_alltoall:', myrank, ierr
!      endif
!      deallocate(qTota2as)


      elseif (ireducemethod == 1) then ! file based approach
!
!       Now for each process on M parts push the fields into an i,j pair named
!       ascii file  restartMap.i.j   1<=i<=N  1<=j<=M.  This is embarassingly
!       parallel with no data collision but writes a lot of files
!
        do l=1,irankN
          if(icountN(l).gt.0) then  !  I have data to write for rank l
!
!           generate the filename since this is a file we have to write
!           note this avoids opening files we don't need but each rank does open
!           all the files it will need to write to....lots of files=weakness
!
!
!           write the number of vtx being mapped from part j to part i
!           note ndof will need to change to the total size of vector field being mapped
!           and this also should be computed in some routine above for allocation
!           purposes.  currently it is at risk of being confused with ndof
!
            itmp=1
            if (l .gt. 0) itmp = int(log10(float(l)))+1
            write (fmt1,"('(''restartMap.'',i',i1,',1x)')") itmp
            write (fnamer,fmt1) l
            fnamer = trim(fnamer) // cname2(myrank+1)
            if (isbinary == 1) then
              open(unit=700+l,file=fnamer,status='unknown',
     &                              form='unformatted')
            else
              open(unit=700+l,file=fnamer,status='unknown')
            endif
          endif
        enddo
!
!       so that we can read directly into the mapped array later, it is helpful to
!       find the number of verts on each part in the N partition which is easily
!       obtained  by looking at the max of each vtx id going to each N rank
!       and then doing an allReduce to find global maximum
!
!       Write solution
!
        do l=1,irankN
          if(icountN(l).gt.0) then  !  I have data to write for rank l
            if (isbinary == 1) then
               write(700+l) icountN(l), ndof
            else
               write(700+l,*) icountN(l), ndof
            endif
          endif
        enddo

        do i=1,nshg
          ifileN=pID(i)+1.1    !shift pID to start from 1
          intvID=(vID(i)+1.1)  !shift vID to start from 1
          if (isbinary == 1) then
            write(700+ifileN) intvID,(qold(i,j),j=1,ndof)
          else
            write(700+ifileN,925) intvID,(qold(i,j),j=1,ndof)
          endif
!         the above "dealing" of records to files based on pID is safe because each
!         process 1..M has a unique set of files to write to
        enddo

!
!       Write dwal
!
        ione = 1
        do l=1,irankN
          if(icountN(l).gt.0) then  !  I have data to write for rank l
            if (isbinary == 1) then
               write(700+l) icountN(l), ione
            else
               write(700+l,*) icountN(l), ione
            endif
          endif
        enddo

        do i=1,nshg
          ifileN=pID(i)+1.1    !shift pID to start from 1
          intvID=(vID(i)+1.1)  !shift vID to start from 1
          if (isbinary == 1) then
            write(700+ifileN) intvID, dwal(i)
          else
            write(700+ifileN,925) intvID, dwal(i)
          endif
!         the above "dealing" of records to files based on pID is safe because each
!         process 1..M has a unique set of files to write to
        enddo
!
!       close all the files on each of the M ranks and deallocate as we are done
!       with the ascii file creation step
!
        do l=1,irankN
          if(icountN(l).gt.0) close(700+l)
        enddo

      endif ! ireducemethod

!
!  Deallocate some memory
!
      call mpi_barrier(mpi_comm_world, ierr)
      if(myrank == 0) then
        write(*,*) 'Deallocating some memory'
      endif
      deallocate(vID)
      deallocate(pID)
      deallocate(qold)
      if(iybar == 1) then
        deallocate(ybar)
      endif
      if(ierror == 1) then
        deallocate(errors)
      endif
      if(numphavg .gt. 0) then
        deallocate(yphbar)
      endif
      if(ivort == 1) then
        deallocate(vort)
      endif
      if(idwal == 1) then
        deallocate(dwal)
      endif
!
!  Make sure all the ranks are done with their restartMap files before
!  trying to read them for ireducemethod = 1
!
      call mpi_barrier(mpi_comm_world, ierr)

!
!  Create a subcommunicator for the first N ranks
!

!      islesseqN = 0
!      if(myrank.lt.irankN) then
!        islesseqN = 1
!      endif
!      call mpi_comm_split(mpi_comm_world, islesseqN, myrank,
!     &                    my_local_comm, ierr)

!      my_local_size = -1
!      call mpi_comm_size(my_local_comm, my_local_size, ierr)
!      my_local_rank = -1
!      call mpi_comm_rank(my_local_comm, my_local_rank, ierr)
!      write(*,*) "RANKS:", myrank, my_local_rank, my_local_size, irankN

!
!     Now for all N ranks, load the partial files written above and merge them
!
!      if(myrank.lt.irankN) then   ! only parts 1..N do work in this phase

        if(ireducemethod == 1) then ! Read the files back for the file based approach.
                                    ! Beware, this method is not maintained any more

          l=myrank+1                ! help me find my mapped files
          do j=1,numpe              ! potentially all M have a map for me
!
!           create the filename for all M of my potential Maps
!
            itmp=1
            if (l .gt. 0) itmp = int(log10(float(l)))+1
            write (fmt1,"('(''restartMap.'',i',i1,',1x)')") itmp
            write (fnamer,fmt1) l
            fnamer = trim(fnamer) // cname2(j)
!
!           most won't exist so inquire for existance
!
            inquire(file=fnamer,exist=exmap)
            if(exmap) then  ! so this map does exist so we open the Map files

              if(isbinary == 1) then
                open(unit=700+l,file=fnamer,status='unknown',
     &                                    form='unformatted')
              else
                open(unit=700+l,file=fnamer,status='unknown')
              endif
!
!             Read solution first
!
              if(isbinary == 1) then
                read(700+l) icountRead, ndofRead  ! read it's header
              else
                read(700+l,*) icountRead, ndofRead  ! read it's header
              endif
              allocate(qreadN(ndofRead))
              do i=1, icountRead                  ! loop over its contents
                if(isbinary == 1) then
                  read(700+l) intvID,(qreadN(k),k=1,ndofRead) ! read a line in binary format
                else
                  read(700+l,925) intvID,(qreadN(k),k=1,ndofRead) ! read a line in ascii format
                endif
! DEBUG
!              if((i.eq.1).and.(myrank.eq.0))
!     &          write(*,*) vidRead, ndofRead, (qreadN(k),k=1,ndofRead)
! DEBUG
                qoldN(intvID,1:ndof)=qreadN(1:ndof) ! fill directly into qoldN
!
!               the above avoids sorting and sifting as duplicated vtx values just overwrite
!               and you get the last one as a final value.  Assuming the field was commu-ed
!               all duplicat values are the same so this is safe.  For error fields this
!               might not be true but it should be o.k.
!
! DEBUG
!              if((i.eq.1).and.(myrank.eq.0))
!     &        write(*,*) intvID,(qoldN(intvID,k),k=1,ndofRead)
! DEBUG
              enddo ! end of this map file so close the file and look for more
              deallocate(qreadN)
!
!             Read dwalN second
!
              if(isbinary == 1) then
                read(700+l) icountRead, ndofRead  ! read it's header
              else
                read(700+l,*) icountRead, ndofRead  ! read it's header
              endif
              do i=1, icountRead
                if(isbinary == 1) then
                  read(700+l) intvID,qreadN1
                else
                  read(700+l,925) intvID,qreadN1
                endif
                dwalN(intvID)=qreadN1 ! fill directly into dwalN
              enddo ! end of this map file so close the file and look for more

              close(700+l)
            endif ! map existed
          enddo ! loop over all M ranks

        endif ! ireducemethod = 1
!
!       At this point I should have a complete field for my part so I write the
!       mapped file
!

! DEBUG
!        itmp=1
!        if (lstep .gt. 0) itmp = int(log10(float(lstep)))+1
!        write (fmt1,"('(''restartMapped.'',i',i1,',1x)')") itmp
!        write (fnamer,fmt1) lstep
!        fnamer = trim(fnamer) // cname2(myrank+1)
!        open(unit=700,file=fnamer,status='unknown') !,
!        write(700,*) maxvIDout(myrank+1),ndof
!        do i=1,maxVIDout(myrank+1)+1
!          write(700,927) i, (qoldN(i,j),j=1,ndof)
!          enddo
!        close(700)
! DEBUG

!
!      Write the N posix files. Note that the nodes on the part boundaries still need to be updated
!
!      Posix format
!        call write_restartonly(myrank, lstep, nshgN, ndof, qoldN)
!        call write_field (myrank,'a','dwal',4,dwalN,'d',nshgN,1,lstep)
!        call write_field (myrank,'a','errors',6,errorsN,'d',
!     &                                        nshgN,ndoferrors,lstep)
!        call write_field (myrank,'a','ybar',4,ybarN,'d',
!     &                                           nshgN,ndofybar,lstep)

!        SyncIO format


!
!      Test the result from mpi_alltoall
!
!      qoldN=-9.87654321e32
!      l=myrank+1
!      do j=1,numpe ! should really go to irankN
!        do i=1,maxicountNglob
!          intvID = qTota2ar( (i-1)*numvar+1 , j) + 1.1
!          if (intvID > 0) then
!            do idof=1,ndof
!              qoldN(intvID,idof) = qTota2ar( (i-1)*numvar+1+idof, j)
!            enddo
!          endif
!        enddo
!      enddo

!      itmp=1
!      if (lstep .gt. 0) itmp = int(log10(float(lstep)))+1
!      write (fmt1,"('(''restartMapped2.'',i',i1,',1x)')") itmp
!      write (fnamer,fmt1) lstep
!      fnamer = trim(fnamer) // cname2(myrank+1)
!      open(unit=700,file=fnamer,status='unknown') !,
!      write(700,*) maxvIDout(myrank+1),ndof
!      do i=1,maxVIDout(myrank+1)+1
!        write(700,927) i, (qoldN(i,j),j=1,ndof)
!      enddo
!      close(700)

!      endif ! my rank is within range 1..N

      call mpi_barrier(mpi_comm_world, ierr)
      if(myrank == 0) then
        write(*,*) 'Writing the reduced restartRedTmp-dat files'
      endif

!      call Write_M2N(myrank, irankN, lstep, nshgN,
!     &          ndof, ndofybar, ndoferrors,
!     &          qoldN, ybarN, errorsN, dwalN)


      nsynciofieldswriterestart = 1+iybar+ierror+numphavg+ivort+idwal

!     We check here how many geombcRed files have been red.
!     If 0, then two possibilities:
!       - Reduction to 1 part so geombcRed are not required. Set nsynciofilesred to 1
!       - User forgot to create the link to geombcRed.
!         Then, assumes the number of geombc and geombcRed files ia the same

      if(nsynciofilesred == 0) then
        if(irankN == 1) then ! 1 part reduction -> set to 1
          nsynciofilesred = 1
        else ! Links are missing
          nsynciofilesred = nsynciofiles
        endif
      endif

      call Write_M2N_SolOnly(myrank,irankN,lstep,nshgN,ndof,qoldN)
      deallocate(qoldN)

      if(iybar == 1) then
        call Write_M2N_Field(myrank,irankN,'a','ybar',4,ybarN,'d',
     &                    nshgN,ndofybar,lstep)
        deallocate(ybarN)
      endif

      if(ierror == 1) then
        call Write_M2N_Field(myrank,irankN,'a','errors',6,errorsN,'d',
     &                    nshgN,ndoferrors,lstep)
        deallocate(errorsN)
      endif

      if(numphavg .gt. 0) then
        do iphavg=1,numphavg
          call write_M2N_phavg2(myrank,irankN,'a','phase_average',13,
     &             iphavg,numphavg,yphbarN(:,:,iphavg),'d',
     &             nshgN,ndofyphbar,lstep)
        enddo
        deallocate(yphbarN)
      endif

      if(ivort == 1) then
        call Write_M2N_Field(myrank,irankN,'a','vorticity',9,
     &                    vortN,'d',nshgN,ndofvort,lstep)
        deallocate(vortN)
      endif

      if(idwal == 1) then
        call Write_M2N_Field(myrank,irankN,'a','dwal',4,dwalN,'d',
     &                    nshgN,1,lstep)
        deallocate(dwalN)
      endif


!
!     Below is debug stuff trying to figure out why the code crashes on exit
!

!       call mpi_barrier(mpi_comm_world, ierr)
!       write(*,*) myrank,'leaving readnblk'

!925   format(i5,2x,6(e14.7,2x))
925   format(i9,1x,6(e24.17,1x)) ! used for file based approach

926   format(9(i5,2x))  !used for debugging
927   format(i5,2x,6(e14.7,2x))  !used for debugging

      return
!
!
      end