ic.f

Go to the documentation of this file.

c-----------------------------------------------------------------------
      subroutine setics
C-----------------------------------------------------------------------
C
C     Set initial conditions.
C
C-----------------------------------------------------------------------
      include 'SIZE'
      include 'DEALIAS'
      include 'INPUT'
      include 'IXYZ'
      include 'GEOM'
      include 'SOLN'
      include 'MASS'
      include 'MVGEOM'
      include 'PARALLEL'
      include 'TSTEP'
 
      logical  iffort(  ldimt1,0:lpert)
     $       , ifrest(0:ldimt1,0:lpert)
     $       , ifprsl(  ldimt1,0:lpert)
 
      LOGICAL  IFANYP
      common /inelr/ nelrr
      common /ctmp1/ work(lx1,ly1,lz1,lelv)
     $ ,             ta1(lx2,ly1,lz1)
     $ ,             ta2(lx2,ly2,lz1)
      integer*8 ntotg,nn
 
      real psmax(ldimt)
 
      if(nio.eq.0) write(6,*) 'set initial conditions'
 
 
      nxyz2=lx2*ly2*lz2       ! Initialize all fields:
      ntot2=nxyz2*nelv
      nxyz1=lx1*ly1*lz1
      ntott=nelt*nxyz1
      ntotv=nelv*nxyz1
      ltott=lelt*nxyz1
 
      call rzero(vx,ntott)
      call rzero(vy,ntott)
      call rzero(vz,ntott)
      call rzero(pr,nxyz2*nelt)
      do 10 ifld=1,ldimt
         call rzero(t(1,1,1,1,ifld),ntott)
   10 continue
 
      jp = 0                  ! Set counter for perturbation analysis
 
      irst = param(46)        ! for lee's restart (rarely used)
      if (irst.gt.0)  call setup_convect(2)
 
c     If moving geometry then add a perturbation to the
c     mesh coordinates (see Subroutine INIGEOM)
 
      if (ifmvbd) call ptbgeom
 
C     Find out what type of i.c. is requested
C     Current options: 
C
C     (1) - User specified fortran function (default is zero i.c.)
C     (2) - Restart from file(s)
C     (3) - Activate pre-solver => steady diffusion / steady Stokes
C
C     If option (2) is requested, also return with the name of the
C     restart file(s) together with the associated dump number
 
      call slogic (iffort,ifrest,ifprsl,nfiles)
 
C      ***** TEMPERATURE AND PASSIVE SCALARS ******
C
C     Check if any pre-solv necessary for temperature/passive scalars
 
      ifanyp = .false.
      DO 100 ifield=2,nfield
         IF (ifprsl(ifield,jp)) THEN
            IF (nio.EQ.0) WRITE(6,101) ifield
            ifanyp = .true.
         ENDIF
  100 CONTINUE
  101 FORMAT(2x,'Using PRESOLVE option for field',i2,'.')
 
C     Fortran function initial conditions for temp/pass. scalars.
      maxfld = nfield
      if (ifmhd) maxfld = npscal+3
 
c     Always call nekuic (pff, 12/7/11)
      do ifield=1,maxfld
         if (nio.eq.0) write(6,*) 'nekuic (1) for ifld ', ifield
         call nekuic
      enddo
 
C     If any pre-solv, do pre-solv for all temperatur/passive scalar fields
      if (ifanyp) call prsolvt
 
      jp = 0 ! jp=0 --> base field, not perturbation field
      do 200 ifield=2,maxfld
         if (iffort(ifield,jp)) then
            if (nio.eq.0) write(6,*) 'call nekuic for ifld ', ifield
            call nekuic
         endif
 200  continue
 
      if (ifpert) then
         ifield=2
         do jp=1,npert
         if (nio.eq.0) write(6,*) 'nekuicP',ifield,jp,iffort(ifield,jp)
            if (iffort(ifield,jp)) call nekuic
         enddo
      endif
      jp = 0
     
 
      call nekgsync()
      call restart(nfiles) !  Check restart files
      call nekgsync()
 
 
C      ***** VELOCITY ******
C     (If restarting for V, we're done,
C     ...else, do pre-solv for fluid if requested.)
 
      ifield = 1
      if (ifprsl(ifield,jp)) call prsolvv
 
 
C     Fortran function initial conditions for velocity.
      ifield = 1
      if (iffort(ifield,jp)) then
         if (nio.eq.0) write(6,*) 'call nekuic for vel  '
         call nekuic
      endif
c
      if (ifpert) then
         ifield=1
         do jp=1,npert
            if (iffort(ifield,jp)) call nekuic
            if (nio.eq.0) write(6,*) 'ic vel pert:',iffort(1,jp),jp
         enddo
      endif
      jp = 0
 
      ntotv = lx1*ly1*lz1*nelv
 
C     Initial mesh velocities
      if (ifmvbd) call opcopy (wx,wy,wz,vx,vy,vz)
      if (ifmvbd.and..not.ifrest(0,jp)) call meshv (2)
 
C     If convection-diffusion of a passive scalar with a fixed velocity field,
C     make sure to fill up lagged arrays since this will not be done in
C     the time-stepping procedure (no flow calculation) (01/18/91 -EMR).
 
      if (.not.ifflow.and.ifheat) then
         itest=0
         DO 400 ifield=2,nfield
            IF (ifadvc(ifield)) itest=1
 400     CONTINUE
         IF (itest.EQ.1) THEN
            nbdmax = 3
            nbdsav = nbdinp
            nbdinp = nbdmax
            DO 500 i=1,nbdmax
               CALL lagvel
 500        CONTINUE
            nbdinp = nbdsav
         ENDIF
      ENDIF
     
C     Ensure that all processors have the same time as node 0.
      if (nid.ne.0) time=0.0
      time=glsum(time,1)
 
      nxyz1=lx1*ly1*lz1
      ntott=nelt*nxyz1
      ntotv=nelv*nxyz1
      nn = nxyz1
      ntotg=nelgv*nn
 
      if (.not.ifdg) then
         ifield = 2
         if (ifflow) ifield = 1
         call rone(work,ntotv)
         ifield = 1
         call dssum  (work,lx1,ly1,lz1)
         call col2   (work,vmult,ntotv)
         rdif  = glsum(work,ntotv)
         rtotg = ntotg
         rdif  = (rdif-rtotg)/rtotg
         if (abs(rdif).gt.1e-14) then
            if(nid.eq.0)write(*,*)'ERROR: dssum test has failed!',rdif
            call exitt
         endif
      endif
 
      call projfld_c0 ! ensure fields are contiguous
 
C print min values
      xxmax = glmin(xm1,ntott)
      yymax = glmin(ym1,ntott)
      zzmax = glmin(zm1,ntott)
 
      vxmax = glmin(vx,ntotv)
      vymax = glmin(vy,ntotv)
      vzmax = glmin(vz,ntotv)
      prmax = glmin(pr,ntot2)
 
      ntot = nxyz1*nelfld(2)
      ttmax = glmin(t ,ntott)
 
      do i=1,ldimt-1
         ntot = nxyz1*nelfld(i+2)
         psmax(i) = glmin(t(1,1,1,1,i+1),ntot)
      enddo
 
      if (nio.eq.0) then
         write(6,19) xxmax,yymax,zzmax
   19    format(' xyz min  ',5g13.5)
      endif
      if (nio.eq.0) then
         write(6,20) vxmax,vymax,vzmax,prmax,ttmax
   20    format(' uvwpt min',5g13.5)
      endif
      if (ldimt-1.gt.0) then
         if (nio.eq.0) write(6,21) (psmax(i),i=1,ldimt-1)
   21    format(' PS min   ',50g13.5)
      endif
 
c print max values
      xxmax = glmax(xm1,ntott)
      yymax = glmax(ym1,ntott)
      zzmax = glmax(zm1,ntott)
 
      vxmax = glmax(vx,ntotv)
      vymax = glmax(vy,ntotv)
      vzmax = glmax(vz,ntotv)
      prmax = glmax(pr,ntot2)
 
      ntot = nxyz1*nelfld(2)
      ttmax = glmax(t ,ntott)
 
      do i=1,ldimt-1
         ntot = nxyz1*nelfld(i+2)
         psmax(i) = glmax(t(1,1,1,1,i+1),ntot)
      enddo
 
      if (nio.eq.0) then
         write(6,16) xxmax,yymax,zzmax
   16    format(' xyz max  ',5g13.5)
      endif
 
      if (nio.eq.0) then
         write(6,17) vxmax,vymax,vzmax,prmax,ttmax
   17    format(' uvwpt max',5g13.5)
      endif
 
      if (ldimt-1.gt.0) then
         if (nio.eq.0)  then
            write(6,18) (psmax(i),i=1,ldimt-1)
   18       format(' PS max   ',50g13.5)
         endif
      endif
 
      if (iflomach .and. ifdp0dt) then
        if (p0th.le.0) call exitti('Invalid thermodynamic pressure!$',1)
        if (gamma0.lt.0) call exitti('Invalid gamma0!$',1)
      endif
 
      if (ifrest(0,jp)) then !  mesh has been read in.
         if (nio.eq.0) write(6,*) 'Restart: recompute geom. factors.'
         call geom_reset(1)  !  recompute geometric factors
      endif
 
      if(nio.eq.0) then
        write(6,*) 'done :: set initial conditions'
        write(6,*) ' '
      endif
 
      return
      end
c-----------------------------------------------------------------------
      subroutine slogic (iffort,ifrest,ifprsl,nfiles)
C---------------------------------------------------------------------
C
C     Set up logicals for initial conditions.
C
C---------------------------------------------------------------------
      include 'SIZE'
      include 'INPUT'
      include 'RESTART'
c
      logical  iffort(  ldimt1,0:lpert)
     $       , ifrest(0:ldimt1,0:lpert)
     $       , ifprsl(  ldimt1,0:lpert)
c
      character*132 line,fname,cdum
      character*2  s2
      character*1  line1(132)
      equivalence(line1,line)
C
C     Default is user specified fortran function (=0 if not specified)
C
      nfldt = nfield
      if (ifmhd) nfldt = nfield+1
 
      do jp=0,npert
         ifrest(0,jp) = .false.
         do ifld=1,nfldt
            iffort(ifld,jp) = .true.
            ifrest(ifld,jp) = .false.
            ifprsl(ifld,jp) = .false.
         enddo
      enddo
 
      jp = 0
      nfiles=0
 
c     Check for Presolve options     
 
      DO 1000 iline=1,15 
         line=initc(iline)
         CALL capit(line,132)
         IF (indx1(line,'PRESOLV',7).NE.0) THEN
C           found a presolve request
            CALL blank(initc(iline),132)
            CALL ljust(line)
            CALL csplit(cdum,line,' ',1)
C
            IF (ltrunc(line,132).EQ.0) THEN
               IF (nio.EQ.0) WRITE(6,700)
  700          FORMAT(/,2x,'Presolve options: ALL')
C              default - all fields are presolved.
               DO 800 ifield=1,nfldt
                  ifprsl(ifield,jp) = .true.
                  iffort(ifield,jp) = .false.
  800          CONTINUE
            ELSE
C           check line for arguments
C
               ll=ltrunc(line,132)
               IF (nio.EQ.0) WRITE(6,810) (line1(l),l=1,ll)
  810          FORMAT(/,2x,'Presolve options: ',132a1)
C
               IF (indx_cut(line,'U',1).NE.0) THEN
                  ifprsl(1,jp) = .true.
                  iffort(1,jp) = .false.
               ENDIF
C
               IF (indx_cut(line,'T',1).NE.0) THEN
                  ifprsl(2,jp) = .true.
                  iffort(2,jp) = .false.
               ENDIF
C
               DO 900 ifield=3,npscal+2
                  ip=ifield-2
                  WRITE(s2,901) ip
                  IF (indx_cut(line,s2,2).NE.0) THEN
                     ifprsl(ifield,jp) = .true.
                     iffort(ifield,jp) = .false.
                  ENDIF
  900          CONTINUE
  901          FORMAT('P',i1)
            ENDIF
         ENDIF
 1000    CONTINUE
C
C     Check for restart options
C
      jp = 0
      DO 2000 iline=1,15
         if (ifpert) jp=iline-1
         line=initc(iline)
         IF (ltrunc(line,132).NE.0) THEN
C           found a filename
            nfiles=nfiles+1
            initc(nfiles)=line
C
C            IF (NIO.EQ.0.AND.NFILES.EQ.1) WRITE(6,1010) LINE
            IF (nio.EQ.0.) WRITE(6,1010) line
 1010       FORMAT(1x,'Checking restart options: ',a132)
c            IF (NID.EQ.0) WRITE(6,'(A132)') LINE
C
C           Parse restart options
 
            call sioflag(ndumps,fname,line)
 
            if (ifgetx) then
               ifrest(0,jp) = .true.
            endif
            if (ifgetu) then
               iffort(1,jp) = .false.
               ifprsl(1,jp) = .false.
               ifrest(1,jp) = .true.
            endif
            if (ifgett) then
               iffort(2,jp) = .false.
               ifprsl(2,jp) = .false.
               ifrest(2,jp) = .true.
            endif
            do 1900 ifield=3,nfldt
c              write(6,*) 'ifgetps:',(ifgtps(k),k=1,ldimt-1)
               if (ifgtps(ifield-2)) then
                  iffort(ifield,jp) = .false.
                  ifprsl(ifield,jp) = .false.
                  ifrest(ifield,jp) = .true.
               endif
 1900       continue
         endif
 2000 continue
 
      return
      end
c-----------------------------------------------------------------------
      subroutine restart(nfiles)
C----------------------------------------------------------------------
C
C     (1) Open restart file(s)
C     (2) Check previous spatial discretization 
C     (3) Map (K1,N1) => (K2,N2) if necessary
C
C     nfiles > 1 has several implications:
C
C     i.   For std. run, data is taken from last file in list, unless
C          explicitly specified in argument list of filename
C
C     ii.  For MHD and perturbation cases, 1st file is for U,P,T;
C          subsequent files are for B-field or perturbation fields
C
C
C----------------------------------------------------------------------
      include 'SIZE'
      include 'TOTAL'
      include 'RESTART'
 
      common /inelr/ nelrr
 
      parameter(lxr=lx1+6)
      parameter(lyr=ly1+6)
      parameter(lzr=lz1+6)
      parameter(lxyzr=lxr*lyr*lzr)
      parameter(lxyzt=lx1*ly1*lz1*lelt)
      parameter(lpsc9=ldimt+9)
 
      common /scrcg/ pm1(lx1*ly1*lz1,lelv)
      COMMON /scrns/ sdump(lxyzt,7)
      integer mesg(40)
 
C     note, this usage of CTMP1 will be less than elsewhere if NELT ~> 9.
      COMMON /ctmp1/ tdump(lxyzr,lpsc9)
      real*4         tdump
c
      REAL SDMP2(LXYZT,LDIMT)
 
c     cdump comes in via PARALLEL (->TOTAL)
 
      character*30 excoder
      character*1  excoder1(30)
      equivalence(excoder,excoder1)
 
 
      character*132 fname
      character*1  fname1(132)
      equivalence(fname1,fname)
 
      integer       hnami (30)
      character*132 hname
      character*1   hname1(132)
      equivalence(hname,hname1)
      equivalence(hname,hnami )
 
      CHARACTER*132 header
 
C     Local logical flags to determine whether to copy data or not.
      logical ifok,iffmat
      integer iposx,iposz,iposu,iposw,iposp,ipost,ipsps(ldimt1)
 
      logical ifbytsw, if_byte_swap_test
      real*4   bytetest
 
c
      ifok=.false.
      ifbytsw = .false.
 
      if(nfiles.lt.1) return
 
      if(nio.eq.0) write(6,*) 'Reading checkpoint data '
 
c use new reader (only binary support)
      p67 = abs(param(67))
      if (p67.eq.6.0) then
         do ifile=1,nfiles
            call sioflag(ndumps,fname,initc(ifile))
            call mfi(fname,ifile)
         enddo
         call bcast(time,wdsize)! Sync time across processors
         return
      endif
 
c use old reader (for ASCII + old binary support)
      
      if (param(67).lt.1.0) then  ! zero only. should be abs.
         iffmat=.true.  ! ascii
      else
         iffmat=.false. ! binary
      endif
 
      do 6000 ifile=1,nfiles
        call sioflag(ndumps,fname,initc(ifile))
        ierr = 0
        if (nid.eq.0) then
 
          if (iffmat) then
            open (unit=91,file=fname,status='old',err=500)
          else
            len= ltrunc(fname,79)
            call izero (hnami,20)
            call chcopy(hname1,fname,len)
c           test for presence of file
            open (unit=91,file=hname
     $           ,form='unformatted',status='old',err=500)
            close(unit=91)
            call byte_open(hname,ierr)
            if(ierr.ne.0) goto 500
          ENDIF
          ifok = .true.
        endif
 
  500   continue
        call lbcast(ifok)
        if (.not.ifok) goto 5000
         
         ndumps = 1
C
C        Only NODE 0 reads from the disk.
C
         DO 1000 idump=1,ndumps
 
            IF (nid.EQ.0) THEN
                ! read header
               if (iffmat) then
                 ierr = 2
                 if(mod(param(67),1.0).eq.0) then ! old header format
                   if(nelgt.lt.10000) then
                     read(91,91,err=10,end=10)
     $              neltr,nxr,nyr,nzr,rstime,istepr,(excoder1(i),i=1,30)
   91                format(4i4,1x,g13.4,i5,1x,30a1)
                     ierr=0
                   else
                     read(91,92,err=10,end=10)
     $              neltr,nxr,nyr,nzr,rstime,istepr,(excoder1(i),i=1,30)
   92                format(i10,3i4,1p1e18.9,i9,1x,30a1)
                     ierr=0
                   endif
                 else                          ! new head format
                   read(91,'(A132)',err=10,end=10) header
                   read(header,*)
     &                  neltr,nxr,nyr,nzr,rstime,istepr,excoder
                   ierr=0
                 endif
               else
                 if(mod(param(67),1.0).eq.0) then  ! old header format
                   call byte_read(hnami,20,ierr)
                   if(ierr.ne.0) goto 10
                   icase = 2
                   if (nelgt.lt.10000) icase = 1
                   ipass = 0
   93              continue  ! test each possible case  UGLY (7/31/07)
                   if(ipass.lt.2) then
                     ipass = ipass+1
                     if(icase.eq.1) then
                       read(hname,'(4i4,1x,g13.4,i5,1x,30a1)',
     $                     err=94,end=94) 
     $                     neltr,nxr,nyr,nzr,rstime,istepr,
     $                     (excoder1(i),i=1,30)
                       goto 95
                     else
                       read(hname,'(i10,3i4,1P1e18.9,i9,1x,30a1)',
     $                     err=94,end=94)
     $                     neltr,nxr,nyr,nzr,rstime,istepr,
     $                     (excoder1(i),i=1,30)
                       goto 95
                     endif
   94                icase = 3-icase  !  toggle: 2-->1  1-->2
                     goto 93
                   else
                     ierr=2
                     goto 10
                   endif
   95              continue
                 else                         ! new head format
                   call byte_read(header,20,ierr)
                   if(ierr.ne.0) goto 10
                   read(header,*)
     &             neltr,nxr,nyr,nzr,rstime,istepr,excoder
                 endif
                 call byte_read(bytetest,1,ierr)
c                call byte_read2(bytetest,1,ierr)
                 if(ierr.ne.0) goto 10
                 ifbytsw = if_byte_swap_test(bytetest,ierr)
                 if(ierr.ne.0) goto 10
               endif
               mesg(1) = neltr
               mesg(2) = nxr
               mesg(3) = nyr
               mesg(4) = nzr
               write(6,*)  'Read mode: ', param(67)
               write(6,333)'neltr,nxr,nyr,nzr: ', neltr,nxr,nyr,nzr
  333          format(a,i9,3i4)
               call chcopy(mesg(5),excoder1,20)
               len  = 14*isize
            endif
   10       call err_chk(ierr,'Error reading restart header. Abort.$')
 
            IF (idump.EQ.1) THEN
               len  = 14*isize
               call bcast(mesg,len)
               neltr = mesg(1)
               nxr   = mesg(2)
               nyr   = mesg(3)
               nzr   = mesg(4)
               call   chcopy(excoder1,mesg(5),20)
 
               call lbcast(ifbytsw)
 
C              Bounds checking on mapped data.
               IF (nxr.GT.lxr) THEN
                  WRITE(6,20) nxr,lx1
   20             FORMAT(//,2x,
     $            'ABORT:  Attempt to map from',i3,
     $            ' to N=',i3,'.',/,2x,
     $            'NEK5000 currently supports mapping from N+6 or less.'
     $            ,/,2x,'Increase N or LXR in IC.FOR.')
                  CALL exitt
               ENDIF
C
C              Figure out position of data in file "IFILE"
C
               nouts=0
               iposx=0
               iposy=0
               iposz=0
               iposu=0
               iposv=0
               iposw=0
               iposp=0
               ipost=0
               DO 40 i=1,npscal
                  ipsps(i)=0
   40          CONTINUE
 
               ips = 0
               nps = 0
               DO 50 i=1, 30
                  IF (excoder1(i).EQ.'X') THEN
                     nouts=nouts + 1
                     iposx=nouts
                     nouts=nouts+1
                     iposy=nouts
                     IF (if3d) THEN
                        nouts=nouts + 1
                        iposz=nouts
                     ENDIF
                  ENDIF
                  IF (excoder1(i).EQ.'U') THEN
                     nouts=nouts + 1
                     iposu=nouts
                     nouts=nouts+1
                     iposv=nouts
                     IF (if3d) THEN
                        nouts=nouts + 1
                        iposw=nouts
                     ENDIF
                  ENDIF
                  IF (excoder1(i).EQ.'P') THEN
                     nouts=nouts + 1
                     iposp=nouts
                  ENDIF
                  IF (excoder1(i).EQ.'T') THEN
                     nouts=nouts + 1
                     ipost=nouts
                  ENDIF
                  IF(mod(param(67),1.0).eq.0.0) THEN
                    i1 = i1_from_char(excoder1(i))
                    if (0.lt.i1.and.i1.lt.10) then
                       if (i1.le.ldimt1) then
                          nouts=nouts + 1
                          ipsps(i1)=nouts
                       else
                          if (nid.eq.0) write(6,2) i1,i,excoder1(i)
   2                      format(2i4,a1,' PROBLEM W/ RESTART DATA')
                       endif
                    endif
                  ELSE
                    IF(excoder1(i).EQ.'S') THEN
                       READ(excoder1(i+1),'(I1)') nps1
                       READ(excoder1(i+2),'(I1)') nps0
                       nps=10*nps1 + nps0 
                       DO is = 1, nps
                         nouts=nouts + 1
                         ipsps(is)=nouts
                       ENDDO
                       GOTO 50
                    ENDIF
                  ENDIF
   50          CONTINUE
 
               IF (nps.GT.(ldimt-1)) THEN
                  IF (nid.EQ.0) THEN 
                    WRITE(*,'(A)') 
     &               'ERROR: restart file has a NSPCAL > LDIMT'
                    WRITE(*,'(A,I2)') 
     &               'Change LDIMT in SIZE'
                  ENDIF
                  CALL exitt
               ENDIF
 
               lname=ltrunc(fname,132)
               if (nio.eq.0) write(6,61) (fname1(i),i=1,lname)
               if (nio.eq.0) write(6,62) 
     $             iposu,iposv,iposw,iposp,ipost,nps,nouts
   61          FORMAT(/,2x,'Restarting from file ',132a1)
   62          FORMAT(2x,'Columns for restart data U,V,W,P,T,S,N: ',7i4)
 
C              Make sure the requested data is present in this file....
               if (iposx.eq.0) ifgetx=.false.
               if (iposy.eq.0) ifgetx=.false.
               if (iposz.eq.0) ifgetz=.false.
               if (iposu.eq.0) ifgetu=.false.
               if (iposv.eq.0) ifgetu=.false.
               if (iposw.eq.0) ifgetw=.false.
               if (iposp.eq.0) ifgetp=.false.
               if (ipost.eq.0) ifgett=.false.
               do 65 i=1,npscal
                  if (ipsps(i).eq.0) ifgtps(i)=.false.
   65          continue
 
C              End of restart file header evaluation.
 
            ENDIF
C
C           Read the error estimators
C           not supported at the moment => just do dummy reading
C
            ifok = .false.
            IF(nid.EQ.0)THEN
               if (iffmat)
     &             READ(91,'(6G11.4)',END=15)(CDUMP,I=1,NELTR)
               ifok = .true.
            ENDIF
  15        continue
            call lbcast(ifok)
            if(.not.ifok) goto 1600
C
C           Read the current dump, double buffer so that we can
C           fit the data on a distributed memory machine,
C           and so we won't have to read the restart file twice
C           in case of an incomplete data file.
C
            nxyzr = nxr*nyr*nzr
C
C           Read the data
C
            nelrr = min(nelgt,neltr) ! # of elements to _really_read
                                     ! why not just neltr?
            do 200 ieg=1,nelrr
               ifok = .false.
               IF (nid.EQ.0) THEN
                 IF (mod(ieg,10000).EQ.1) WRITE(6,*) 'Reading',ieg
                 IF (iffmat) THEN
                    READ(91,*,err=70,END=70)
     $              ((tdump(ixyz,ii),ii=1,nouts),ixyz=1,nxyzr)
                 ELSE
                    do ii=1,nouts
                       call byte_read(tdump(1,ii),nxyzr,ierr)
                       if(ierr.ne.0) then
                          write(6,*) "Error reading xyz restart data"
                          goto 70
                       endif
                    enddo
                 ENDIF
                 ifok=.true.
               ENDIF
C
C              Notify other processors that we've read the data OK.
C
  70           continue
               call lbcast(ifok)
               IF (.NOT.ifok) GOTO 1600
C
C              MAPDMP maps data from NXR to lx1
C              (and sends data to the appropriate processor.)
C
C              The buffer SDUMP is used so that if an incomplete dump
C              file is found (e.g. due to UNIX io buffering!), then
C              the previous read data stored in VX,VY,.., is not corrupted.
C
               IF (ifgetx) CALL mapdmp
     $         (sdump(1,1),tdump(1,iposx),ieg,nxr,nyr,nzr,ifbytsw)
               IF (ifgetx) CALL mapdmp
     $         (sdump(1,2),tdump(1,iposy),ieg,nxr,nyr,nzr,ifbytsw)
               IF (ifgetz) CALL mapdmp
     $         (sdump(1,3),tdump(1,iposz),ieg,nxr,nyr,nzr,ifbytsw)
               IF (ifgetu) CALL mapdmp
     $         (sdump(1,4),tdump(1,iposu),ieg,nxr,nyr,nzr,ifbytsw)
               IF (ifgetu) CALL mapdmp
     $         (sdump(1,5),tdump(1,iposv),ieg,nxr,nyr,nzr,ifbytsw)
               IF (ifgetw) CALL mapdmp
     $         (sdump(1,6),tdump(1,iposw),ieg,nxr,nyr,nzr,ifbytsw)
               IF (ifgetp) CALL mapdmp
     $         (sdump(1,7),tdump(1,iposp),ieg,nxr,nyr,nzr,ifbytsw)
               if (ifgett) call mapdmp
     $         (sdmp2(1,1),tdump(1,ipost),ieg,nxr,nyr,nzr,ifbytsw)
 
C              passive scalars
               do 100 ips=1,npscal
                  if (ifgtps(ips)) call mapdmp(sdmp2(1,ips+1)
     $               ,tdump(1,ipsps(ips)),ieg,nxr,nyr,nzr,ifbytsw)
  100          continue
 
  200       continue
C
C           Successfully read a complete field, store it.
C
            nerr = 0              ! Count number of elements rec'd by nid
            do ieg=1,nelrr
               mid = gllnid(ieg)
               if (mid.eq.nid) nerr = nerr+1
            enddo
 
            nxyz2=lx2*ly2*lz2
            nxyz1=lx1*ly1*lz1
            ntott=nerr*nxyz1
            ntotv=nerr*nxyz1   ! Problem for differing Vel. and Temp. counts!
                               ! for now we read nelt dataset
 
            if (ifmhd.and.ifile.eq.2) then
               if (ifgetu) call copy(bx,sdump(1,4),ntott)
               if (ifgetu) call copy(by,sdump(1,5),ntott)
               if (ifgetw) call copy(bz,sdump(1,6),ntott)
               if (ifgetp) then
                 if (nio.eq.0) write(6,*) 'getting restart pressure'
                 if (ifsplit) then
                    call copy( pm,sdump(1,7),ntotv)
                 else
                  do iel=1,nelv
                    iiel = (iel-1)*nxyz1+1
                    call map12 (pm(1,1,1,iel),sdump(iiel,7),iel)
                  enddo
                 endif
               endif
               if (ifaxis.and.ifgett) 
     $            call copy(t(1,1,1,1,2),sdmp2(1,1),ntott)
            elseif (ifpert.and.ifile.ge.2) then
               j=ifile-1  ! pointer to perturbation field
               if (ifgetu) call copy(vxp(1,j),sdump(1,4),ntotv)
               if (ifgetu) call copy(vyp(1,j),sdump(1,5),ntotv)
               if (ifgetw) call copy(vzp(1,j),sdump(1,6),ntotv)
               if (ifgetp) then
                  if (nio.eq.0) write(6,*) 'getting restart pressure'
                  if (ifsplit) then
                     call copy(prp(1,j),sdump(1,7),ntotv)
                  else
                     do ie=1,nelv
                        ie1 = (ie-1)*nxyz1+1
                        ie2 = (ie-1)*nxyz2+1
                        call map12 (prp(ie2,j),sdump(ie1,7),ie)
                     enddo
                  endif
               endif
               if (ifgett) call copy(tp(1,1,j),sdmp2(1,1),ntott)
C              passive scalars
               do ips=1,npscal
                  if (ifgtps(ips))
     $            call copy(tp(1,ips+1,j),sdmp2(1,ips+1),ntott)
               enddo
 
            else  ! Std. Case
               if (ifgetx) call copy(xm1,sdump(1,1),ntott)
               if (ifgetx) call copy(ym1,sdump(1,2),ntott)
               if (ifgetz) call copy(zm1,sdump(1,3),ntott)
               if (ifgetu) call copy(vx ,sdump(1,4),ntotv)
               if (ifgetu) call copy(vy ,sdump(1,5),ntotv)
               if (ifgetw) call copy(vz ,sdump(1,6),ntotv)
               if (ifgetp) call copy(pm1,sdump(1,7),ntotv)
               if (ifgett) call copy(t,sdmp2(1,1),ntott)
C              passive scalars
               do i=1,npscal
                  if (ifgtps(i))
     $            call copy(t(1,1,1,1,i+1),sdmp2(1,i+1),ntott)
               enddo
 
               if (ifaxis) call axis_interp_ic(pm1)      ! Interpolate to axi mesh
               if (ifgetp) call map_pm1_to_pr(pm1,ifile) ! Interpolate pressure
 
               if (ifgtim) time=rstime
            endif
 
 1000    CONTINUE
         GOTO 1600
 
 1600    CONTINUE
C
         IF (idump.EQ.1.AND.nid.EQ.0) THEN
            write(6,1700) fname
            write(6,1701) ieg,ixyz
            write(6,1702) 
     $            ((tdump(jxyz,ii),ii=1,nouts),jxyz=ixyz-1,ixyz)
 1700       FORMAT(5x,'WARNING:  No data read in for file ',a132)
 1701       FORMAT(5x,'Failed on  element',i4,',  point',i5,'.')
 1702       FORMAT(5x,'Last read dump:',/,5g15.7)
            write(6,*) nid,'call exitt 1702a',idump
            call exitt
         ELSEIF (idump.EQ.1) THEN
            write(6,*) nid,'call exitt 1702b',idump
            call exitt
         ELSE
            idump=idump-1
            IF (nio.EQ.0) WRITE(6,1800) idump
 1800       FORMAT(2x,'Successfully read data from dump number',i3,'.')
         ENDIF
         if (iffmat) then
            if (nid.eq.0) close(unit=91)
         else
            ierr = 0
            if (nid.eq.0) call byte_close(ierr)
            call err_chk(ierr,'Error closing restart file in restart$')
         endif
         GOTO 6000
 
C        Can't open file...
 5000    CONTINUE
         if (nid.eq.0) write(6,5001) fname 
 5001    FORMAT(2x,'   *******   ERROR   *******    '
     $       ,/,2x,'   *******   ERROR   *******    '
     $       ,/,2x,'   Could not open restart file:'
     $       ,/,a132
     $      ,//,2x,'Quitting in routine RESTART.')
         CLOSE(unit=91)
         call exitt
 5002    CONTINUE
         IF (nio.EQ.0) WRITE(6,5001) hname 
         call exitt
C
C
C     End of IFILE loop
 6000 CONTINUE
C
      return
      end
C
c-----------------------------------------------------------------------
      subroutine sioflag(ndumps,fname,rsopts)
C
C     Set IO flags according to Restart Options File, RSOPTS
C
      include 'SIZE'
      include 'INPUT'
      include 'RESTART'
      include 'TSTEP'
 
      character*132 rsopts,fname
      character*2  s2
      logical ifgtrl
 
C     Scratch variables..
      logical ifdeft,ifanyc
      CHARACTER*132 RSOPT     ,LINE
      CHARACTER*1  RSOPT1(132),LINE1(132)
      equivalence(rsopt1,rsopt)
      equivalence(line1,line)
C
C     Parse filename
C
C        CSPLIT splits S1 into two parts, delimited by S2.  
C        S1 returns with 2nd part of S1.  CSPLIT returns 1st part.
C
      rsopt=rsopts
      call ljust(rsopt)
      call csplit(fname,rsopt,' ',1)
C     check fname for user supplied extension.
      if (indx1(fname,'.',1).eq.0) then
         len=ltrunc(fname,132)
         len1=len+1
         len4=len+4
         fname(len1:len4)='.fld'
      endif
C
C     Parse restart options
C
C     set default flags
C
      ifgetx=.false.
      ifgetz=.false.
      ifgetu=.false.
      ifgetw=.false.
      ifgetp=.false.
      ifgett=.false.
      do 100 i=1,ldimt-1
         ifgtps(i)=.false.
  100 continue
      ifgtim=.true.
      ndumps=0
C
C     Check for default case - just a filename given, no i/o options specified
C
      ifdeft=.true.
C
C     Parse file for i/o options and/or dump number
C
      CALL capit(rsopt,132)
 
      IF (ltrunc(rsopt,132).NE.0) THEN
C
C        Check for explicit specification of restart TIME.
C
         ito=indx_cut(rsopt,'TIME',4)
         ifgtim=.true.
         IF (ito.NE.0) THEN
C           user has specified the time explicitly.
            it1=indx_cut(rsopt,'=',1)
            it8=132-it1
            CALL blank(line,132)
            CALL chcopy(line,rsopt1(it1),it8)
            IF (ifgtrl(ttime,line)) THEN
               ifgtim=.false.
               time=ttime
            ENDIF
C           remove the user specified time from the RS options line.
            ita=132-ito+1
            CALL blank(rsopt1(ito),ita)
            CALL ljust(line)
            it1=indx1(line,' ',1)
            itb=132-it1+1
            CALL chcopy(rsopt1(ito),line1(it1),itb)
         ENDIF
 
C        Parse field specifications.
 
         ixo=indx_cut(rsopt,'X',1)
         IF (ixo.NE.0) THEN
            ifdeft=.false.
            ifgetx=.true.
            IF (if3d) ifgetz=.true.
         ENDIF
 
         ivo=indx_cut(rsopt,'U',1)
         IF (ivo.NE.0) THEN
            ifdeft=.false.
            ifgetu=.true.
            IF (if3d) ifgetw=.true.
         ENDIF
 
         ipo=indx_cut(rsopt,'P',1)
         IF (ipo.NE.0) THEN
            ifdeft=.false.
            ifgetp=.true.
         ENDIF
 
         ito=indx_cut(rsopt,'T',1)
         IF (ito.NE.0) THEN
            ifdeft=.false.
            ifgett=.true.
         ENDIF
 
         do 300 i=1,ldimt-1
            write (s2,301) i
            ipo=indx_cut(rsopt,s2,2)
            if (ipo.ne.0) then
               ifdeft=.false.
               ifgtps(i)=.true.
            endif
  300    continue
  301    format('S',i1)
 
C        Get number of dumps from remainder of user supplied line.
         if (ifgtrl(tdumps,rsopt)) ndumps=int(tdumps)
      endif
 
C     If no fields were explicitly specified, assume getting all fields. 
      if (ifdeft) then
         ifgetx=.true.
         IF (if3d) ifgetz=.true.
         ifanyc=.false.
         DO 400 i=1,nfield
            IF (ifadvc(i)) ifanyc=.true.
  400    CONTINUE
         IF (ifflow.OR.ifanyc) THEN
            ifgetu=.true.
            IF (if3d) ifgetw=.true.
         ENDIF
         if (ifflow) ifgetp=.true.
         if (ifheat) ifgett=.true.
#ifdef CMTNEK
         ifgett=.true. ! CMT-nek still not compatible with IFHEAT
#endif
         do 410 i=1,ldimt-1
            ifgtps(i)=.true.
  410    continue
      endif
 
 
 
 
      return
      END
c-----------------------------------------------------------------------
      subroutine mapdmp(sdump,tdump,ieg,nxr,nyr,nzr,if_byte_sw)
C----------------------------------------------------------------------
C
C----------------------------------------------------------------------
      include 'SIZE'
      include 'PARALLEL'
C
      parameter(lxyz1=lx1*ly1*lz1)
      parameter(lxr=lx1+6)
      parameter(lyr=ly1+6)
      parameter(lzr=lz1+6)
      parameter(lxyzr=lxr*lyr*lzr)
C
      REAL   SDUMP(LXYZ1,LELT)
      real*4 tdump(lxyzr)
c
      logical if_byte_sw
c
      nxyz=lx1*ly1*lz1
      nxyr=nxr*nyr*nzr
      ierr=0
C
C     Serial processor code:
C
      IF (np.EQ.1) THEN
 
         IF (if_byte_sw) call byte_reverse(tdump,nxyr,ierr)
         if(ierr.ne.0) call exitti("Error in mapdmp")
         IF (nxr.EQ.lx1.AND.nyr.EQ.ly1.AND.nzr.EQ.lz1) THEN
            CALL copy4r(sdump(1,ieg),tdump,nxyz)
         ELSE
C           do the map    (assumes that NX=NY=NZ, or NX=NY, NZ=1)
            call mapab4r(sdump(1,ieg),tdump,nxr,1)
         ENDIF
 
      ELSE
C
C     Parallel code - send data to appropriate processor and map.
C
         jnid=gllnid(ieg)
         mtype=3333+gllel(ieg)
         len=4*nxyr
         le1=4
         IF (nid.EQ.0.AND.jnid.NE.0) THEN
c           hand-shake
            CALL csend(mtype,tdump,le1,jnid,nullpid)
            CALL crecv2(mtype,dummy,le1,jnid)
            CALL csend(mtype,tdump,len,jnid,nullpid)
         ELSEIF (nid.NE.0.AND.jnid.EQ.nid) THEN
C           Receive data from node 0
            CALL crecv2(mtype,dummy,le1,0)
            CALL csend(mtype,tdump,le1,0,nullpid)
            CALL crecv2(mtype,tdump,len,0)
         ENDIF
C
C        If the data is targeted for this processor, then map 
C        to appropriate element.
C
         IF (jnid.EQ.nid) THEN
            ie=gllel(ieg)
            IF (if_byte_sw) call byte_reverse(tdump,nxyr,ierr)
            IF (nxr.EQ.lx1.AND.nyr.EQ.ly1.AND.nzr.EQ.lz1) THEN
               CALL copy4r(sdump(1,ie),tdump,nxyz)
            ELSE
               call mapab4r(sdump(1,ie),tdump,nxr,1)
            ENDIF
         ENDIF
         call err_chk(ierr,'Error using byte_reverse in mapdmp.$')
C
C        End of parallel distribution/map routine.
C
      ENDIF
      return
      END
c-----------------------------------------------------------------------
      subroutine mapab(x,y,nxr,nel)
C---------------------------------------------------------------
C
C     Interpolate Y(NXR,NYR,NZR,NEL) to X(lx1,ly1,lz1,NEL)
C     (assumes that NXR=NYR=NZR, or NXR=NYR, NZR=1)
C---------------------------------------------------------------
C
      include 'SIZE'
      include 'IXYZ'
      include 'WZ'
C
      parameter(lxr=lx1+6)
      parameter(lyr=ly1+6)
      parameter(lzr=lz1+6)
      parameter(lxyzr=lxr*lyr*lzr)
      parameter(lxyz1=lx1*ly1*lz1)
      dimension x(lx1,ly1,lz1,nel)
      dimension y(nxr,nxr,nxr,nel)
 
      common /ctmp0/  xa(lxyzr)      ,xb(lx1,ly1,lzr) ,xc(lxyzr)
     $              , zgmr(lxr)      ,wgtr(lxr)
      common /ctmpab/ ires(lxr,lxr)  ,itres(lxr,lxr)
      real ires,itres
 
      INTEGER NOLD
      SAVE    nold
      DATA    nold /0/
 
C     Bounds checking on mapped data.
      if (nxr.gt.lxr) then
         if (nid.eq.0) write(6,20) nxr,lx1
   20    FORMAT(//,2x,
     $   'ABORT:  Attempt to map from',i3,
     $   ' to N=',i3,'.',/,2x,
     $   'NEK5000 currently supports mapping from N+6 or less.'
     $   ,/,2x,'Increase N')
         call exitt
      endif
 
      nzr = nxr
      IF(lz1.EQ.1) nzr=1
      nyzr = nxr*nzr
      nxy1 = lx1*ly1
 
      IF (nxr.NE.nold) THEN
         nold=nxr
         CALL zwgll   (zgmr,wgtr,nxr)
         CALL igllm   (ires,itres,zgmr,zgm1,nxr,lx1,nxr,lx1)      
         IF (nio.EQ.0) WRITE(6,10) nxr,lx1
   10       FORMAT(2x,'Mapping restart data from Nold=',i2
     $               ,' to Nnew=',i2,'.')
      ENDIF
C
      DO 1000 ie=1,nel
         CALL mxm (ires,lx1,y(1,1,1,ie),nxr,xa,nyzr)
         DO 100 iz=1,nzr
            izoff = 1 + (iz-1)*lx1*nxr
            CALL mxm (xa(izoff),lx1,itres,nxr,xb(1,1,iz),ly1)
  100    CONTINUE
         IF (ldim.EQ.3) THEN
            CALL mxm (xb,nxy1,itres,nzr,x(1,1,1,ie),lz1)
         ELSE
            CALL copy(x(1,1,1,ie),xb,nxy1)
         ENDIF
 1000 CONTINUE
C
      return
      END
c-----------------------------------------------------------------------
      subroutine mapab4r(x,y,nxr,nel)
C---------------------------------------------------------------
C
C     Interpolate Y(NXR,NYR,NZR,NEL) to X(lx1,ly1,lz1,NEL)
C     (assumes that NXR=NYR=NZR, or NXR=NYR, NZR=1)
c
c     Input:  real*4,  Output:  default precision
c
C---------------------------------------------------------------
C
      include 'SIZE'
      include 'IXYZ'
      include 'WZ'
C
      parameter(lxr=lx1+6)
      parameter(lyr=ly1+6)
      parameter(lzr=lz1+6)
      parameter(lxyzr=lxr*lyr*lzr)
      parameter(lxyz1=lx1*ly1*lz1)
      real*4 x(lx1,ly1,lz1,nel)
      REAL   Y(NXR,NXR,NXR,NEL)
 
      common /ctmp0/  xa(lxyzr)      ,xb(lx1,ly1,lzr) ,xc(lxyzr)
     $              , zgmr(lxr)      ,wgtr(lxr)
      common /ctmpa4/ ires(lxr,lxr)  ,itres(lxr,lxr)
      real ires,itres
 
      INTEGER NOLD
      SAVE    nold
      DATA    nold /0/
 
C     Bounds checking on mapped data.
      if (nxr.gt.lxr) then
         if (nid.eq.0) write(6,20) nxr,lx1
   20    FORMAT(//,2x,
     $   'ABORT:  Attempt to map from',i3,
     $   ' to N=',i3,'.',/,2x,
     $   'NEK5000 currently supports mapping from N+6 or less.'
     $   ,/,2x,'Increase N')
         call exitt
      endif
 
      nzr = nxr
      IF(lz1.EQ.1) nzr=1
      nyzr = nxr*nzr
      nxy1 = lx1*ly1
      nxyzr = nxr*nxr*nzr
 
      IF (nxr.NE.nold) THEN
         nold=nxr
         CALL zwgll   (zgmr,wgtr,nxr)
         CALL igllm   (ires,itres,zgmr,zgm1,nxr,lx1,nxr,lx1)      
         IF (nio.EQ.0) WRITE(6,10) nxr,lx1
   10       FORMAT(2x,'Mapping restart data from Nold=',i2
     $               ,' to Nnew=',i2,'.')
      ENDIF
C
      DO 1000 ie=1,nel
         call copy4r(xc,y(1,1,1,ie),nxyzr)
         CALL mxm (ires,lx1,xc,nxr,xa,nyzr)
         DO 100 iz=1,nzr
            izoff = 1 + (iz-1)*lx1*nxr
            CALL mxm (xa(izoff),lx1,itres,nxr,xb(1,1,iz),ly1)
  100    CONTINUE
         IF (ldim.EQ.3) THEN
            CALL mxm (xb,nxy1,itres,nzr,x(1,1,1,ie),lz1)
         ELSE
            CALL copy(x(1,1,1,ie),xb,nxy1)
         ENDIF
 1000 CONTINUE
C
      return
      END
c-----------------------------------------------------------------------
      function i1_from_char(s1)
      character*1 s1
 
      character*10 n10
      save         n10
      data         n10 / '0123456789' /
 
      i1_from_char = indx2(n10,10,s1,1)-1
 
      return
      end
c-----------------------------------------------------------------------
      integer function indx2(s1,l1,s2,l2)
      character*132 s1,s2
 
      n1=l1-l2+1
      indx2=0
      if (n1.lt.1) return
 
      do i=1,n1
         i2=i+l2-1
         if (s1(i:i2).eq.s2(1:l2)) then
            indx2=i
            return
         endif
      enddo
 
      return
      end
c-----------------------------------------------------------------------
      INTEGER FUNCTION indx1(S1,S2,L2)
      CHARACTER*132 s1,s2
C
      n1=132-l2+1
      indx1=0
      IF (n1.LT.1) return
C
      DO 100 i=1,n1
         i2=i+l2-1
         IF (s1(i:i2).EQ.s2(1:l2)) THEN
            indx1=i
            return
         ENDIF
  100 CONTINUE
C
      return
      END
c-----------------------------------------------------------------------
      INTEGER FUNCTION indx_cut(S1,S2,L2)
C
C     INDX_CUT is returned with the location of S2 in S1 (0 if not found)
C     S1     is returned with 1st occurance of S2 removed.
C
      CHARACTER*1 s1(132),s2(132)
C
      i1=indx1(s1,s2,l2)
C
      IF (i1.NE.0) THEN
C
         n1=132-l2
         DO 100 i=i1,n1
            i2=i+l2
C           remove the 1st occurance of S2 from S1.
            s1(i)=s1(i2)
  100    CONTINUE
         n2=n1+1
         DO 200 i=n2,132
            s1(i)=' '
  200    CONTINUE
      ENDIF
C
      indx_cut=i1
      return
      END
c-----------------------------------------------------------------------
      subroutine csplit(s0,s1,s2,l0)
      CHARACTER*132 S0,S1,S2
C     split string S1 into two parts, delimited by S2.
C
      i2=indx_cut(s1,s2,l0)
      IF (i2.EQ.0) return
C
      i1=i2-1
      CALL blank(s0,132)
      s0(1:i1)=s1(1:i1)
      CALL lshft(s1,i2)
C
      return
      END
c-----------------------------------------------------------------------
      subroutine lshft(string,ipt)
C     shift string from IPT to the left
C     INPUT : "abcde......    test    "
C     OUTPUT: "e......    test        "     if ipt.eq.5
      CHARACTER*1 STRING(132)
C
      DO 20 j=1,133-ipt
         ij=ipt+j-1
         string(j)=string(ij)
   20 CONTINUE
      DO 30 j=134-ipt,132
         string(j)=' '
   30 CONTINUE
      return
      END
c-----------------------------------------------------------------------
      subroutine ljust(string)
C     left justify string
      CHARACTER*1 STRING(132)
C
      IF (string(1).NE.' ') return
C
      DO 100 i=2,132
C
         IF (string(i).NE.' ') THEN
            DO 20 j=1,133-i
               ij=i+j-1
               string(j)=string(ij)
   20       CONTINUE
            DO 30 j=134-i,132
               string(j)=' '
   30       CONTINUE
            return
         ENDIF
C
  100 CONTINUE
      return
      END
c-----------------------------------------------------------------------
      subroutine chknorm (ifzero)
C----------------------------------------------------------------------
C
C     Check if trivial user specified initial conditions
C
C----------------------------------------------------------------------
      include 'SIZE'
      include 'INPUT'
      include 'TSTEP'
      LOGICAL IFZERO
C
      ifzero = .true.
C
      IF (ifflow) THEN
         ifield = 1
         imesh  = 1
         CALL unorm
         IF (vnrml8.GT.0.) ifzero = .false.
      ENDIF
      IF (ifheat) THEN
         DO 100 ifield=2,nfield
            imesh = 1
            IF (iftmsh(ifield)) imesh = 2
            CALL unorm
            IF (tnrml8(ifield).GT.0.) ifzero = .false.
 100     CONTINUE
      ENDIF
c
      return
      END
C
c-----------------------------------------------------------------------
      subroutine prsolvt
C----------------------------------------------------------------------
C
C     Use steady state solution as initial condition 
C     for temperatur/passive scalar
C
C----------------------------------------------------------------------
      include 'SIZE'
      include 'INPUT'
      include 'TSTEP'
      LOGICAL  IFSAV1,IFSAV2(LDIMT1)
C
      IF (nio.EQ.0) WRITE(6,*) ' '
      IF (nio.EQ.0) WRITE(6,*) 'Conduction pre-solver activated'
C
C     Set logical IFTRAN to false (steady state)
C     Save logicals for convection
C     Turn convection off
C
      ifsav1 = iftran
      iftran = .false.
      DO 100 ifield=2,nfield
         ifsav2(ifield) = ifadvc(ifield)
         ifadvc(ifield) = .false.
 100  CONTINUE
C
      CALL setprop
      CALL setsolv
C
      IF(nio.EQ.0)WRITE(6,*)'Steady conduction/passive scalar problem'
C
      DO 200 igeom=1,2
         CALL heat (igeom)
 200  CONTINUE
C
C     Set IFTRAN to true again
C     Turn convection on again
C
      iftran = ifsav1
      DO 300 ifield=2,nfield
         ifadvc(ifield) = ifsav2(ifield)
 300  CONTINUE
C
      return
      END
C
c-----------------------------------------------------------------------
      subroutine prsolvv
C----------------------------------------------------------------------
C
C     Use steady Stokes solution as initial condition 
C     for flow problem
C
C----------------------------------------------------------------------
      include 'SIZE'
      include 'INPUT'
      include 'SOLN'
      include 'TSTEP'
      LOGICAL  IFSAV1,IFSAV2
C
      IF (nio.EQ.0) WRITE(6,*) ' '
      IF (nio.EQ.0) WRITE(6,*) 'Velocity pre-solver activated'
C
C     Initialize velocity to some non-trivial RHS to avoid FP trap in i860.
C
      IF (param(60).NE.0.0) THEN
         small=10.0e-10
         CALL cfill(vx,small,ntotv)
         CALL cfill(vy,small,ntotv)
         CALL cfill(vz,small,ntotv)
      ENDIF
C
C     Set logical IFTRAN to false (steady state)
C     Save logicals for convection
C     Turn convection off
C
      IF (ifsplit) THEN
        WRITE(6,10)
   10   FORMAT(
     $ /,2x,'ERROR: Steady Stokes Flow initial condition cannot'
     $,/,2x,'       be computed using the splitting formulation.'
     $,/,2x,'       Either compute using UZAWA, or remove PRESOLVE'
     $,/,2x,'       request for velocity.'
     $,/,2x
     $,/,2x,'       ABORTING IN PRSOLVV.')
        CALL exitt
      ENDIF
C
      ifsav1 = iftran
      ifsav2 = ifadvc(ifield)
      iftran = .false.
      ifadvc(ifield) = .false.
C
      CALL setprop
      CALL setsolv
C
      IF (nio.EQ.0) WRITE (6,*) 'Steady Stokes problem'
      DO 100 igeom=1,2
         IF (.NOT.ifsplit) CALL fluid (igeom)
 100  CONTINUE
C
C     Set IFTRAN to true again
C     Turn convection on again
C
      iftran = ifsav1
      ifadvc(ifield) = ifsav2
C
      return
      END
C
c-----------------------------------------------------------------------
      subroutine nekuic ! user specified fortran function (=0 if not specified)
 
      include 'SIZE'
      include 'INPUT'
      include 'SOLN'
      include 'TSTEP'
      include 'PARALLEL'
      include 'NEKUSE'
 
      integer e,eg
 
      nel   = nelfld(ifield)
 
      do e=1,nel
         eg = lglel(e)
         do 300 k=1,lz1
         do 300 j=1,ly1
         do 300 i=1,lx1
           call nekasgn (i,j,k,e)
           call useric  (i,j,k,eg)
           if (jp.eq.0) then
             if (ifield.eq.1) then
               vx(i,j,k,e) = ux
               vy(i,j,k,e) = uy
               vz(i,j,k,e) = uz
             elseif (ifield.eq.ifldmhd .and. ifmhd) then
               bx(i,j,k,e) = ux
               by(i,j,k,e) = uy
               bz(i,j,k,e) = uz
             else
               t(i,j,k,e,ifield-1) = temp
             endif
           else
             ijke = i+lx1*((j-1)+ly1*((k-1) + lz1*(e-1)))
             if (ifield.eq.1) then
               vxp(ijke,jp) = ux
               vyp(ijke,jp) = uy
               vzp(ijke,jp) = uz
             else
               tp(ijke,ifield-1,jp) = temp
             endif
           endif
 
 300     continue
      enddo
 
      return
      END
c-----------------------------------------------------------------------
      subroutine capit(lettrs,n)
C     Capitalizes string of length n
      CHARACTER LETTRS(N)
C
      DO 5 i=1,n
         int=ichar(lettrs(i))
         IF(int.GE.97 .AND. int.LE.122) THEN
            int=int-32
            lettrs(i)=char(int)
         ENDIF
5     CONTINUE
      return
      END
c-----------------------------------------------------------------------
      LOGICAL FUNCTION ifgtrl(VALUE,LINE)
C
C     Read VALUE from LINE and set IFGTRL to .TRUE. if successful,
C                                  IFGTRL to .FALSE. otherwise.
C
C     This complicated function is necessary thanks to the Ardent,
C     which won't allow free formatted reads (*) from internal strings!
C
      CHARACTER*132 line
      CHARACTER*132 work
      CHARACTER*8  fmat
C
C     Note that the format Fn.0 is appropriate for fields of type:
C          34   34.0  34.0e+00
C     The only difficulty would be with '34' but since we identify
C     the field width exactly, there is no problem.
C
      ifgtrl=.false.
      VALUE=0.0
C
      work=line
      CALL ljust(work)
      ifldw=indx1(work,' ',1)-1
C
      IF (ifldw.GT.0) THEN
         WRITE(fmat,10) ifldw
   10    FORMAT('(F',i3.3,'.0)')
         READ(work,fmat,err=100,END=100) tval
         VALUE=tval
         ifgtrl=.true.
         return
      ENDIF
C
  100 CONTINUE
      return
      END
c-----------------------------------------------------------------------
      LOGICAL FUNCTION ifgtil(IVALUE,LINE)
C
C     Read IVALUE from LINE and set IFGTIL to .TRUE. if successful,
C                                   IFGTIL to .FALSE. otherwise.
C
C     This complicated function is necessary thanks to the Ardent,
C     which won't allow free formatted reads (*) from internal strings!
C
      CHARACTER*132 line
      CHARACTER*132 work
      CHARACTER*8  fmat
C
      ifgtil=.false.
      ivalue=0
C
      work=line
      CALL ljust(work)
      ifldw=indx1(work,' ',1)-1
C
      IF (ifldw.GT.0) THEN
         WRITE(fmat,10) ifldw
   10    FORMAT('(F',i3.3,'.0)')
         READ(work,fmat,err=100,END=100) tval
         ivalue=int(tval)
         ifgtil=.true.
         return
      ENDIF
C
  100 CONTINUE
      return
      END
c-----------------------------------------------------------------------
      subroutine perturb(tt,ifld,eps)
      include 'SIZE'
      include 'TOTAL'
c
      real tt(1)
      integer ifld
 
      ifield = ifld
 
      n = lx1*ly1*lz1*nelfld(ifield)
      call vcospf(tt,bm1,n)
      call cmult(tt,eps,n)
      call dssum(tt,lx1,ly1,lz1)
 
      return
      end
c-----------------------------------------------------------------------
      subroutine vcospf(x,y,n)
      real x(1),y(1)
      do i=1,n
         x(i) = cos(1000.*y(i))
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine vbyte_swap(x,n)
      character*1 x(0:3,1),tmp0,tmp1
      character*1 in (0:3), out(0:3)
      real*4      in4     , out4
      equivalence(in ,in4 )
      equivalence(out,out4)
c
      do i=1,n
         do j=0,3
            in(j) = x(j,i)
         enddo
         tmp0   = x(0,i)
         tmp1   = x(1,i)
         x(0,i) = x(3,i)
         x(1,i) = x(2,i)
         x(2,i) = tmp1
         x(3,i) = tmp0
         do j=0,3
            out(j) = x(j,i)
         enddo
         write(6,*) 'swap:',i,in4,out4
      enddo
c
      return
      end
c-----------------------------------------------------------------------
      logical function if_byte_swap_test(bytetest,ierr)
      include 'SIZE'
 
      real*4 bytetest,test2
      real*4 test_pattern
      save   test_pattern
 
      test_pattern = 6.54321
      eps          = 0.00020
      etest        = abs(test_pattern-bytetest)
      if_byte_swap_test = .true.
      if (etest.le.eps) if_byte_swap_test = .false.
 
      test2 = bytetest
      call byte_reverse(test2,1,ierr)
      if (nid.eq.0 .and. loglevel.gt.2) 
     $   write(6,*) 'byte swap:',if_byte_swap_test,bytetest,test2
      return
      end
c-----------------------------------------------------------------------
      subroutine geom_reset(icall)
C
C     Generate geometry data
C
      include 'SIZE'
      include 'INPUT'
      include 'GEOM'
      include 'SOLN'
      include 'TSTEP'
      include 'WZ'
c
      COMMON /scruz/ xm3(lx1,ly1,lz1,lelt)
     $ ,             ym3(lx1,ly1,lz1,lelt)
     $ ,             zm3(lx1,ly1,lz1,lelt)
C
c
      if(nio.eq.0) write(6,*) 'regenerate geometry data',icall
 
      ntot = lx1*ly1*lz1*nelt
c
      if (lx3.eq.lx1) then
         call copy(xm3,xm1,ntot)
         call copy(ym3,ym1,ntot)
         call copy(zm3,zm1,ntot)
      else
         call map13_all(xm3,xm1)
         call map13_all(ym3,ym1)
         if (if3d) call map13_all(zm3,zm1)
      endif
c
      CALL geom1 (xm3,ym3,zm3)
      CALL geom2
      CALL updmsys (1)
      CALL volume
      CALL setinvm
      CALL setdef
      CALL sfastax
c
      do ie = 1,nelt
         ! x
         xc(1,ie) = xm1(1  ,1  ,1  ,ie)
         xc(2,ie) = xm1(lx1,1  ,1  ,ie)
         xc(3,ie) = xm1(lx1,ly1,1  ,ie)
         xc(4,ie) = xm1(1  ,ly1,1  ,ie)
         xc(5,ie) = xm1(1  ,1  ,lz1,ie)
         xc(6,ie) = xm1(lx1,1  ,lz1,ie)
         xc(7,ie) = xm1(lx1,ly1,lz1,ie)
         xc(8,ie) = xm1(1  ,ly1,lz1,ie)
         ! y
         yc(1,ie) = ym1(1  ,1  ,1  ,ie)
         yc(2,ie) = ym1(lx1,1  ,1  ,ie)
         yc(3,ie) = ym1(lx1,ly1,1  ,ie)
         yc(4,ie) = ym1(1  ,ly1,1  ,ie)
         yc(5,ie) = ym1(1  ,1  ,lz1,ie)
         yc(6,ie) = ym1(lx1,1  ,lz1,ie)
         yc(7,ie) = ym1(lx1,ly1,lz1,ie)
         yc(8,ie) = ym1(1  ,ly1,lz1,ie)
         ! z 
         zc(1,ie) = zm1(1  ,1  ,1  ,ie)
         zc(2,ie) = zm1(lx1,1  ,1  ,ie)
         zc(3,ie) = zm1(lx1,ly1,1  ,ie)
         zc(4,ie) = zm1(1  ,ly1,1  ,ie)
         zc(5,ie) = zm1(1  ,1  ,lz1,ie)
         zc(6,ie) = zm1(lx1,1  ,lz1,ie)
         zc(7,ie) = zm1(lx1,ly1,lz1,ie)
         zc(8,ie) = zm1(1  ,ly1,lz1,ie)
      enddo
 
      if(nio.eq.0) then
        write(6,*) 'done :: regenerate geometry data',icall
        write(6,*) ' '
      endif
 
      return
      end
c-----------------------------------------------------------------------
      subroutine dsavg(u)
c
c
      include 'SIZE'
      include 'TOTAL'
      real u(lx1,ly1,lz1,lelt)
c
c
c     Take direct stiffness avg of u
c
c
      ifieldo = ifield
      if (ifflow) then
         ifield = 1
         ntot = lx1*ly1*lz1*nelv
         call dssum(u,lx1,ly1,lz1)
         call col2 (u,vmult,ntot)
      else
         ifield = 2
         ntot = lx1*ly1*lz1*nelt
         call dssum(u,lx1,ly1,lz1)
         call col2 (u,tmult,ntot)
      endif
      ifield = ifieldo
c
      return
      end
c-----------------------------------------------------------------------
      subroutine map13_all(x3,x1)
c
      include 'SIZE'
      include 'TOTAL'
c
      real x3(lx3,ly3,lz3,lelt)
      real x1(lx1,ly1,lz1,lelt)
c
      integer e
c
      do e=1,nelt
         call map13 (x3(1,1,1,e),x1(1,1,1,e),e)
      enddo
c
      return
      end
c-----------------------------------------------------------------------
      subroutine mfi_gets(u,wk,lwk,iskip)
 
      include 'SIZE'
      include 'INPUT'
      include 'PARALLEL'
      include 'RESTART'
 
 
      real u(lx1*ly1*lz1,1)
 
      real*4 wk(2*lwk) ! message buffer
      parameter(lrbs=20*lx1*ly1*lz1*lelt)
      common /vrthov/ w2(lrbs) ! read buffer
      real*4 w2
 
      integer e,ei,eg,msg_id(lelt)
      logical iskip
      integer*8 i8tmp
 
      call nekgsync() ! clear outstanding message queues.
 
      nxyzr  = nxr*nyr*nzr  
      dnxyzr = nxyzr 
      len    = nxyzr*wdsizr  ! message length
      if (wdsizr.eq.8) nxyzr = 2*nxyzr
 
      ! check message buffer wk
      num_recv  = nxyzr*nelt
      num_avail = size(wk)
      call lim_chk(num_recv,num_avail,'     ','     ','mfi_gets a')
 
      ! setup read buffer
      if (nid.eq.pid0r) then
         i8tmp = int(nxyzr,8)*int(nelr,8)
         nread = i8tmp/int(lrbs,8)
         if (mod(i8tmp,int(lrbs,8)).ne.0) nread = nread + 1
         if(ifmpiio) nread = iglmax(nread,1) ! needed because of collective read
         nelrr = nelr/nread
      endif
      call bcast(nelrr,4)
      call lim_chk(nxyzr*nelrr,lrbs,'     ','     ','mfi_gets b')
 
      ! pre-post recieves
      if (np.gt.1) then
         l = 1
         do e=1,nelt
            msg_id(e) = irecv(e,wk(l),len)
            l = l+nxyzr
         enddo
      endif
 
      ierr = 0
      if (nid.eq.pid0r.and.np.gt.1) then ! only i/o nodes will read
         ! read blocks of size nelrr
         k = 0
         do i = 1,nread
            if(i.eq.nread) then ! clean-up 
              nelrr = nelr - (nread-1)*nelrr
              if(nelrr.lt.0) nelrr = 0
            endif
            
            if(ierr.eq.0) then
              if(ifmpiio) then
                call byte_read_mpi(w2,nxyzr*nelrr,-1,ifh_mbyte,ierr)
              else
                call byte_read (w2,nxyzr*nelrr,ierr)
              endif
            endif
 
            ! distribute data across target processors
            l = 1
            do e = k+1,k+nelrr
               jnid = gllnid(er(e))                ! where is er(e) now?
               jeln = gllel(er(e))
               if(ierr.ne.0) call rzero(w2(l),len)
               call csend(jeln,w2(l),len,jnid,0)  ! blocking send
               l = l+nxyzr
            enddo
            k  = k + nelrr
         enddo
      elseif (np.eq.1) then
         if(ifmpiio) then
           call byte_read_mpi(wk,nxyzr*nelr,-1,ifh_mbyte,ierr)
         else
           call byte_read(wk,nxyzr*nelr,ierr)
         endif
      endif
 
      if (iskip) then
         call nekgsync() ! clear outstanding message queues.
         goto 100     ! don't use the data
      endif
 
      nxyzr = nxr*nyr*nzr
      nxyzv = nxr*nyr*nzr
      nxyzw = nxr*nyr*nzr
      if (wdsizr.eq.8) nxyzw = 2*nxyzw
 
      l = 1
      do e=1,nelt
         if (np.gt.1) then
            call msgwait(msg_id(e))
            ei = e
         elseif(np.eq.1) then
            ei = er(e)
         endif
         if (if_byte_sw) then
            if(wdsizr.eq.8) then
              call byte_reverse8(wk(l),nxyzv*2,ierr)
            else
              call byte_reverse(wk(l),nxyzv,ierr)
            endif
         endif
         if (nxr.eq.lx1.and.nyr.eq.ly1.and.nzr.eq.lz1) then
            if (wdsizr.eq.4) then         ! COPY
               call copy4r(u(1,ei),wk(l        ),nxyzr)
            else
               call copy  (u(1,ei),wk(l        ),nxyzr)
            endif
         else                             ! INTERPOLATE
            if (wdsizr.eq.4) then
               call mapab4r(u(1,ei),wk(l        ),nxr,1)
            else
               call mapab  (u(1,ei),wk(l        ),nxr,1)
            endif
         endif
         l = l+nxyzw
      enddo
 
 
 100  call err_chk(ierr,'Error reading restart data,in gets.$')
      return
      end
c-----------------------------------------------------------------------
      subroutine mfi_getv(u,v,w,wk,lwk,iskip)
 
      include 'SIZE'
      include 'INPUT'
      include 'PARALLEL'
      include 'RESTART'
 
      real u(lx1*ly1*lz1,1),v(lx1*ly1*lz1,1),w(lx1*ly1*lz1,1)
      logical iskip
 
      real*4 wk(2*lwk) ! message buffer
      parameter(lrbs=20*lx1*ly1*lz1*lelt)
      common /vrthov/ w2(lrbs) ! read buffer
      real*4 w2
 
      integer e,ei,eg,msg_id(lelt)
      integer*8 i8tmp
 
      call nekgsync() ! clear outstanding message queues.
 
      nxyzr  = ldim*nxr*nyr*nzr
      len    = nxyzr*wdsizr             ! message length in bytes
      if (wdsizr.eq.8) nxyzr = 2*nxyzr
 
      ! check message buffer wk
      num_recv  = nxyzr*nelt 
      num_avail = size(wk)
      call lim_chk(num_recv,num_avail,'     ','     ','mfi_getv a')
 
      ! setup read buffer
      if(nid.eq.pid0r) then
         i8tmp = int(nxyzr,8)*int(nelr,8)
         nread = i8tmp/int(lrbs,8)
         if (mod(i8tmp,int(lrbs,8)).ne.0) nread = nread + 1
         if(ifmpiio) nread = iglmax(nread,1) ! needed because of collective read
         nelrr = nelr/nread
      endif
      call bcast(nelrr,4)
      call lim_chk(nxyzr*nelrr,lrbs,'     ','     ','mfi_getv b')
 
      ! pre-post recieves (one mesg per element)
      ! this assumes we never pre post more messages than supported
      if (np.gt.1) then
         l = 1
         do e=1,nelt
            msg_id(e) = irecv(e,wk(l),len)
            l = l+nxyzr
         enddo
      endif
 
      ierr = 0
      if (nid.eq.pid0r .and. np.gt.1) then ! only i/o nodes
         k = 0
         do i = 1,nread
            if(i.eq.nread) then ! clean-up 
              nelrr = nelr - (nread-1)*nelrr
              if(nelrr.lt.0) nelrr = 0
            endif
 
            if(ierr.eq.0) then
              if(ifmpiio) then 
                call byte_read_mpi(w2,nxyzr*nelrr,-1,ifh_mbyte,ierr)
              else
                call byte_read (w2,nxyzr*nelrr,ierr)
              endif
            endif
 
            ! redistribute data based on the current el-proc map
            l = 1
            do e = k+1,k+nelrr
               jnid = gllnid(er(e))                ! where is er(e) now?
               jeln = gllel(er(e))
               if(ierr.ne.0) call rzero(w2(l),len)
               call csend(jeln,w2(l),len,jnid,0)  ! blocking send
               l = l+nxyzr
            enddo
            k  = k + nelrr
         enddo
      elseif (np.eq.1) then
         if(ifmpiio) then 
           call byte_read_mpi(wk,nxyzr*nelr,-1,ifh_mbyte,ierr)
         else
           call byte_read(wk,nxyzr*nelr,ierr)
         endif
      endif
 
      if (iskip) then
         call nekgsync() ! clear outstanding message queues.
         goto 100     ! don't assign the data we just read
      endif
 
      nxyzr = nxr*nyr*nzr
      nxyzv = ldim*nxr*nyr*nzr
      nxyzw = nxr*nyr*nzr
      if (wdsizr.eq.8) nxyzw = 2*nxyzw
 
      l = 1
      do e=1,nelt
         if (np.gt.1) then
            call msgwait(msg_id(e))
            ei = e
         else if(np.eq.1) then
            ei = er(e) 
         endif
         if (if_byte_sw) then
            if(wdsizr.eq.8) then
               call byte_reverse8(wk(l),nxyzv*2,ierr)
            else
               call byte_reverse(wk(l),nxyzv,ierr)
            endif
         endif
         if (nxr.eq.lx1.and.nyr.eq.ly1.and.nzr.eq.lz1) then
            if (wdsizr.eq.4) then         ! COPY
               call copy4r(u(1,ei),wk(l        ),nxyzr)
               call copy4r(v(1,ei),wk(l+  nxyzw),nxyzr)
               if (if3d) 
     $         call copy4r(w(1,ei),wk(l+2*nxyzw),nxyzr)
            else
               call copy  (u(1,ei),wk(l        ),nxyzr)
               call copy  (v(1,ei),wk(l+  nxyzw),nxyzr)
               if (if3d) 
     $         call copy  (w(1,ei),wk(l+2*nxyzw),nxyzr)
            endif
         else                             ! INTERPOLATE
            if (wdsizr.eq.4) then
               call mapab4r(u(1,ei),wk(l        ),nxr,1)
               call mapab4r(v(1,ei),wk(l+  nxyzw),nxr,1)
               if (if3d) 
     $         call mapab4r(w(1,ei),wk(l+2*nxyzw),nxr,1)
            else
               call mapab  (u(1,ei),wk(l        ),nxr,1)
               call mapab  (v(1,ei),wk(l+  nxyzw),nxr,1)
               if (if3d) 
     $         call mapab  (w(1,ei),wk(l+2*nxyzw),nxr,1)
            endif
         endif
         l = l+ldim*nxyzw
      enddo
 
 100  call err_chk(ierr,'Error reading restart data, in getv.$')
      return
      end
c-----------------------------------------------------------------------
      subroutine mfi_parse_hdr(hdr,ierr)
      include 'SIZE'
 
      character*132 hdr
 
      if (indx2(hdr,132,'#std',4).eq.1) then
          call parse_std_hdr(hdr)
      else
         if (nio.eq.0) write(6,80) hdr
         if (nio.eq.0) write(6,80) 'NONSTD HDR, parse_hdr, abort.'
  80     format(a132)
         ierr = 1
      endif
 
      return
      end
c-----------------------------------------------------------------------
      subroutine parse_std_hdr(hdr)
      include 'SIZE'
      include 'INPUT'
      include 'SOLN'
      include 'PARALLEL'
      include 'RESTART'
      include 'TSTEP'
 
      character*132 hdr
      character*4 dummy
      logical if_press_mesh
 
      p0thr = -1
      if_press_mesh = .false.
 
      read(hdr,*,iostat=ierr) dummy
     $         ,  wdsizr,nxr,nyr,nzr,nelr,nelgr,timer,istpr
     $         ,  ifiler,nfiler
     $         ,  rdcode      ! 74+20=94
     $         ,  p0thr, if_press_mesh
 
      if (ierr.gt.0) then ! try again without pressure format flag
        read(hdr,*,iostat=ierr) dummy
     $         ,  wdsizr,nxr,nyr,nzr,nelr,nelgr,timer,istpr
     $         ,  ifiler,nfiler
     $         ,  rdcode      ! 74+20=94
     $         ,  p0thr
      endif
 
      if (ierr.gt.0) then ! try again without mean pressure
        read(hdr,*,err=99) dummy
     $         ,  wdsizr,nxr,nyr,nzr,nelr,nelgr,timer,istpr
     $         ,  ifiler,nfiler
     $         ,  rdcode      ! 74+20=94
      endif
 
c     set if_full_pres flag
      if_full_pres = .false.
      if (.not.ifsplit) if_full_pres = if_press_mesh
 
c      ifgtim  = .true.  ! always get time
      ifgetxr = .false.
      ifgetur = .false.
      ifgetpr = .false.
      ifgettr = .false.
      do k=1,ldimt-1
         ifgtpsr(k) = .false.
      enddo
 
      npsr = 0
      nps  = 0
      do i=1,10 
         if (rdcode1(i).eq.'X') ifgetxr = .true.
         if (rdcode1(i).eq.'U') ifgetur = .true.
         if (rdcode1(i).eq.'P') ifgetpr = .true.
         if (rdcode1(i).eq.'T') ifgettr = .true.
         if (rdcode1(i).eq.'S') then
            read(rdcode1(i+1),'(I1)') nps1
            read(rdcode1(i+2),'(I1)') nps0
            npsr = 10*nps1+nps0
            nps  = npsr
            if(npsr.gt.ldimt-1) nps=ldimt-1
            do k=1,nps
               ifgtpsr(k) = .true.
            enddo
            ! nothing will follow
            GOTO 50
         endif
      enddo
 
  50  if (nps.lt.npsr) then
         if (nid.eq.0) then 
           write(*,'(A,/,A)') 
     &      'WARNING: restart file has a NSPCAL > LDIMT',
     &      'read only part of the fld-data!'
         endif
      endif
 
      if (nps.lt.npscal) then
         if (nid.eq.0) then 
           write(*,'(A,/,A)') 
     &      'WARNING: NPSCAL read from restart file differs from ',
     &      'currently used NPSCAL!'
         endif
      endif
 
      p0th = 1 
      if (p0thr.gt.0) p0th = p0thr
 
      return
 
   99 continue   !  If we got here, then the May 2008 variant of std hdr
                 !  failed and we may have an older input file.
 
      call parse_std_hdr_2006(hdr,rdcode)  ! try the original header format
 
      return
      end
c-----------------------------------------------------------------------
      subroutine parse_std_hdr_2006(hdr,rlcode)
      include 'SIZE'
      include 'INPUT'
      include 'RESTART'
 
      character*132 hdr
      character*1 rlcode(20)
 
c                4  7  10  13   23    33    53    62     68     74
      read(hdr,1) wdsizr,nxr,nyr,nzr,nelr,nelgr,timer,istpr
     $         , ifiler,nfiler
     $         , (rlcode(k),k=1,20)                   ! 74+20=94
    1 format(4x,i2,3i3,2i10,e20.13,i9,2i6,20a1)
 
      if (nid.eq.0) write(6,*) 'WARNING: reading depreacted header!'
 
c     Assign read conditions, according to rdcode
c     NOTE: In the old hdr format: what you see in file is what you get.
c      ifgtim  = .true.  ! always get time
      ifgetxr = .false.
      ifgetur = .false.
      ifgetpr = .false.
      ifgettr = .false.
      do k=1,npscal
         ifgtpsr(k) = .false.
      enddo
 
      if (rlcode(1).eq.'X') ifgetxr = .true.
      if (rlcode(2).eq.'U') ifgetur = .true.
      if (rlcode(3).eq.'P') ifgetpr = .true.
      if (rlcode(4).eq.'T') ifgettr = .true.
      do k=1,npscal
         if (rlcode(4+k).ne.' ') ifgtpsr(k) = .true.
      enddo
 
 
      return
      end
c-----------------------------------------------------------------------
      subroutine mfi(fname_in,ifile)
c
c     (1) Open restart file(s)
c     (2) Check previous spatial discretization 
c     (3) Map (K1,N1) => (K2,N2) if necessary
c
c     nfiles > 1 has several implications:
c
c     i.   For std. run, data is taken from last file in list, unless
c          explicitly specified in argument list of filename
c
c     ii.  For MHD and perturbation cases, 1st file is for U,P,T;
c          subsequent files are for B-field or perturbation fields
c
c
      include 'SIZE'
      include 'TOTAL'
      include 'RESTART'
      character*132  hdr
      character*132  fname_in
 
      character*132  fname
      character*1    fnam1(132)
      equivalence(fnam1,fname)
 
      character*1    frontc
 
      parameter(lwk = 7*lx1*ly1*lz1*lelt)
      common /scrns/ wk(lwk)
      common /scrcg/ pm1(lx1*ly1*lz1,lelv)
      integer e
 
      integer*8 offs0,offs,nbyte,stride,strideB,nxyzr8
 
      tiostart=dnekclock()
 
      ! add full path if required
      call blank(fname,132)
      call chcopy(frontc, fname_in, 1)
      if (frontc .ne. '/') then
        lenp = ltrunc(path,132)
        lenf = ltrunc(fname_in,132)
        call chcopy(fnam1(1),path,lenp)
        call chcopy(fnam1(lenp+1),fname_in,lenf)
      else
        lenf = ltrunc(fname_in,132)
        call chcopy(fnam1(1),fname_in,lenf)     
      endif
 
      call mfi_prepare(fname)       ! determine reader nodes +
                                    ! read hdr + element mapping 
 
      offs0   = iheadersize + 4 + isize*nelgr
      nxyzr8  = nxr*nyr*nzr
      strideb = nelbr* nxyzr8*wdsizr
      stride  = nelgr* nxyzr8*wdsizr
 
      iofldsr = 0
      if (ifgetxr) then      ! if available
         offs = offs0 + ldim*strideb
         call byte_set_view(offs,ifh_mbyte)
         if (ifgetx) then
c            if(nid.eq.0) write(6,*) 'Reading mesh'
            call mfi_getv(xm1,ym1,zm1,wk,lwk,.false.)
         else                ! skip the data
            call mfi_getv(xm1,ym1,zm1,wk,lwk,.true.)
         endif
         iofldsr = iofldsr + ldim
      endif
 
      if (ifgetur) then
         offs = offs0 + iofldsr*stride + ldim*strideb
         call byte_set_view(offs,ifh_mbyte)
         if (ifgetu) then
            if (ifmhd.and.ifile.eq.2) then
c               if(nid.eq.0) write(6,*) 'Reading B field'
               call mfi_getv(bx,by,bz,wk,lwk,.false.)
            else
c               if(nid.eq.0) write(6,*) 'Reading velocity field'
               call mfi_getv(vx,vy,vz,wk,lwk,.false.)
            endif
         else
            call mfi_getv(vx,vy,vz,wk,lwk,.true.)
         endif
         iofldsr = iofldsr + ldim
      endif
 
      if (ifgetpr) then
         offs = offs0 + iofldsr*stride + strideb
         call byte_set_view(offs,ifh_mbyte)
         if (ifgetp) then
c           if(nid.eq.0) write(6,*) 'Reading pressure field'
            call mfi_gets(pm1,wk,lwk,.false.)
         else
            call mfi_gets(pm1,wk,lwk,.true.)
         endif
         iofldsr = iofldsr + 1
      endif
 
      if (ifgettr) then
         offs = offs0 + iofldsr*stride + strideb
         call byte_set_view(offs,ifh_mbyte)
         if (ifgett) then
c            if(nid.eq.0) write(6,*) 'Reading temperature field'
            call mfi_gets(t,wk,lwk,.false.)
         else
            call mfi_gets(t,wk,lwk,.true.)
         endif
         iofldsr = iofldsr + 1
      endif
 
      ierr = 0
      do k=1,ldimt-1
         if (ifgtpsr(k)) then
            offs = offs0 + iofldsr*stride + strideb
            call byte_set_view(offs,ifh_mbyte)
            if (ifgtps(k)) then
c               if(nid.eq.0) write(6,'(A,I2,A)') ' Reading ps',k,' field'
               call mfi_gets(t(1,1,1,1,k+1),wk,lwk,.false.)
            else
               call mfi_gets(t(1,1,1,1,k+1),wk,lwk,.true.)
            endif
            iofldsr = iofldsr + 1
         endif
      enddo
      nbyte = 0
      if(nid.eq.pid0r) nbyte = iofldsr*nelr*wdsizr*nxr*nyr*nzr
 
      if (ifgtim) time = timer
 
      if(ifmpiio) then
        if(nid.eq.pid0r) call byte_close_mpi(ifh_mbyte,ierr)
      else
        if(nid.eq.pid0r) call byte_close(ierr)
      endif
      call err_chk(ierr,'Error closing restart file, in mfi.$')
      tio = dnekclock()-tiostart
 
      dnbyte = nbyte
      nbyte = glsum(dnbyte,1)
      nbyte = nbyte + iheadersize + 4 + isize*nelgr
 
      if (tio.eq.0) tio=1
      if (nio.eq.0) write(6,7) istep,time,
     &             nbyte/tio/1024/1024/10,
     &             nfiler
    7 format(/,i9,1pe12.4,' done :: Read checkpoint data',/,
     &       30x,'avg data-throughput = ',f7.1,'MBps',/,
     &       30x,'io-nodes = ',i5,/)
 
 
      if (ifaxis) call axis_interp_ic(pm1)      ! Interpolate to axi mesh
      if (ifgetp) call map_pm1_to_pr(pm1,ifile) ! Interpolate pressure
 
      return
      end
c-----------------------------------------------------------------------
      subroutine addfid(fname,fid)
      include 'SIZE'
      include 'TSTEP'
      include 'INPUT'
      include 'RESTART'
 
 
      character*1 fname(132)
      integer fid
 
      character*8  eight,fmt,s8
      save         eight
      data         eight / "????????" /
 
      do ipass=1,2      ! 2nd pass, in case 1 file/directory
         do k=8,1,-1
            i1 = indx1(fname,eight,k)
            if (i1.ne.0) then ! found k??? string
               write(fmt,1) k,k
    1          format('(i',i1,'.',i1,')')
               write(s8,fmt) fid
               call chcopy(fname(i1),s8,k)
               goto 10
            endif
         enddo
   10    continue
      enddo
      
      return
      end
c-----------------------------------------------------------------------
      subroutine mfi_prepare(hname)  ! determine which nodes are readers
      character*132 hname
 
      include 'SIZE'
      include 'PARALLEL'
      include 'RESTART'
      include 'INPUT'
 
      integer stride
      character*132 hdr, hname_
      logical if_byte_swap_test
      real*4 bytetest
 
      integer*8 offs0,offs
 
      ierr = 0
      ! rank0 (i/o master) will do a pre-read to get some infos 
      ! we need to have in advance
      if (nid.eq.0) then
         call chcopy(hname_,hname,132)
         call addfid(hname_,0)
         call byte_open(hname_,ierr)
 
         if(ierr.ne.0) goto 101
         call blank     (hdr,iheadersize)
         call byte_read (hdr,iheadersize/4,ierr)
         if(ierr.ne.0) goto 101
         call byte_read (bytetest,1,ierr)
         if(ierr.ne.0) goto 101
         if_byte_sw = if_byte_swap_test(bytetest,ierr) ! determine endianess
         if(ierr.ne.0) goto 101
         call byte_close(ierr)
      endif
 
 101  continue
      call err_chk(ierr,'Error reading restart header in mfi_prepare$')
 
      call bcast(if_byte_sw,lsize) 
      call bcast(hdr,iheadersize)  
      call mfi_parse_hdr(hdr,ierr)
 
      ifmpiio = .false.
      if (nfiler.eq.1 .and. abs(param(67)).eq.6) ifmpiio = .true.
#ifdef NOMPIIO
      ifmpiio = .false.
#endif
 
      if(.not.ifmpiio) then
 
        stride = np / nfiler
        if (stride.lt.1) then
           write(6,*) nfiler,np,'  TOO MANY FILES, mfi_prepare'
           call exitt
        endif
 
        if (mod(nid,stride).eq.0) then ! i/o clients
           pid0r = nid
           pid1r = nid + stride
           fid0r = nid / stride
           call blank(hdr,iheadersize)
 
           call addfid(hname,fid0r)
           if(nid.eq.pid0r) write(6,*) '      FILE:',hname
           call byte_open(hname,ierr)
 
           if(ierr.ne.0) goto 102
           call byte_read (hdr, iheadersize/4,ierr)  
           if(ierr.ne.0) goto 102
           call byte_read (bytetest,1,ierr) 
           if(ierr.ne.0) goto 102
           call mfi_parse_hdr (hdr,ierr)    ! replace hdr with correct one 
           if (nelr.gt.lelr) then
              write(6,*) 'ERROR: increase lelr in SIZE to ', nelr 
              call exitt
           endif
           call byte_read (er,nelr,ierr)    ! get element mapping
           if(if_byte_sw) call byte_reverse(er,nelr,ierr)
        else
           pid0r = 0
           pid1r = 0
           fid0r = 0
        endif
 
      else
 
        pid0r  = nid
        pid1r  = nid
        offs0  = iheadersize + 4
        nfiler = np
 
        ! number of elements to read 
        nelr = nelgr/np
        do i = 0,mod(nelgr,np)-1
           if(i.eq.nid) nelr = nelr + 1
        enddo
        nelbr = igl_running_sum(nelr) - nelr 
        offs = offs0 + nelbr*isize
 
        call addfid(hname,fid0r)
        if(nio.eq.0) write(6,*) '      FILE:',hname
        call byte_open_mpi(hname,ifh_mbyte,.true.,ierr)
 
        if(ierr.ne.0) goto 102
        call byte_set_view(offs,ifh_mbyte)
        if (nelr.gt.lelr) then
           write(6,*) 'ERROR: increase lelr in SIZE to ', nelr 
           call exitt
        endif
        call byte_read_mpi(er,nelr,-1,ifh_mbyte,ierr)
        if(ierr.ne.0) goto 102
        if(if_byte_sw) call byte_reverse(er,nelr,ierr)
 
      endif
 
 102  continue
      call err_chk(ierr,'Error reading header/element map.$')
 
      return
      end
c-----------------------------------------------------------------------
      subroutine axis_interp_ic(pm1)
 
      include 'SIZE'
      include 'TOTAL'
      include 'RESTART'
 
      real pm1(lx1,ly1,lz1,lelv)
 
      common /ctmp0/ axism1(lx1,ly1)
      integer e
 
      if (.not.ifaxis) return
 
      do e=1,nelv
         if (ifrzer(e)) then
           if (ifgetx) then
             call mxm   (xm1(1,1,1,e),lx1,iatlj1,ly1,axism1,ly1)
             call copy  (xm1(1,1,1,e),axism1,lx1*ly1)
             call mxm   (ym1(1,1,1,e),lx1,iatlj1,ly1,axism1,ly1)
             call copy  (ym1(1,1,1,e),axism1,lx1*ly1)
           endif
           if (ifgetu) then
             call mxm    (vx(1,1,1,e),lx1,iatlj1,ly1,axism1,ly1)
             call copy   (vx(1,1,1,e),axism1,lx1*ly1)
             call mxm    (vy(1,1,1,e),lx1,iatlj1,ly1,axism1,ly1)
             call copy   (vy(1,1,1,e),axism1,lx1*ly1)
           endif
           if (ifgetw) then
             call mxm    (vz(1,1,1,e),lx1,iatlj1,ly1,axism1,ly1)
             call copy   (vz(1,1,1,e),axism1,lx1*ly1)
           endif
           if (ifgetp) then
             call mxm    (pm1(1,1,1,e),lx1,iatlj1,ly1,axism1,ly1)
             call copy   (pm1(1,1,1,e),axism1,lx1*ly1)
           endif
           if (ifgett) then
             call mxm  (t(1,1,1,e,1),lx1,iatlj1,ly1,axism1,ly1)
             call copy (t(1,1,1,e,1),axism1,lx1*ly1)
           endif
           do ips=1,npscal
            is1 = ips + 1
            if (ifgtps(ips)) then
             call mxm (t(1,1,1,e,is1),lx1,iatlj1,ly1,axism1,ly1)
             call copy(t(1,1,1,e,is1),axism1,lx1*ly1)
            endif
           enddo
         endif
      enddo
   
      return
      end
c-----------------------------------------------------------------------
      subroutine map_pm1_to_pr(pm1,ifile)
 
      include 'SIZE'
      include 'TOTAL'
      include 'RESTART'
 
      real pm1(lx1*ly1*lz1,lelv)
      integer e
 
      nxyz2 = lx2*ly2*lz2
 
      if (ifmhd.and.ifile.eq.2) then
         do e=1,nelv
            if (if_full_pres) then
               call copy  (pm(1,1,1,e),pm1(1,e),nxyz2)
            else
               call map12 (pm(1,1,1,e),pm1(1,e),e)
            endif
         enddo
      elseif (ifsplit) then
         call copy (pr,pm1,lx1*ly1*lz1*nelv)
      else
         do e=1,nelv
            if (if_full_pres) then
               call copy  (pr(1,1,1,e),pm1(1,e),nxyz2)
            else
               call map12 (pr(1,1,1,e),pm1(1,e),e)
            endif
         enddo
      endif
   
      return
      end
c-----------------------------------------------------------------------
      subroutine full_restart(fnames,n_restart)
      include 'SIZE'
      include 'TOTAL'
 
      character *(*) fnames(*)
 
      ifile = istep+1  ! istep=0,1,...
 
      if (ifile.le.n_restart) then
         p67 = param(67)
         param(67) = 6.00
         call chcopy (initc,fnames(ifile),80)
         call bcast  (initc,80)
         call restart(1)
         call setprop
         param(67)=p67
      endif
   
      return
      end
c-----------------------------------------------------------------------
      subroutine projfld_c0()
 
      include 'SIZE'
      include 'TOTAL'
 
      nxyz1 = lx1*ly1*lz1
      ntott = nelt*nxyz1
      ntotv = nelv*nxyz1
 
      if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'projfld_c0'
 
c     if (ifflow.and..not.ifdg)  then  ! Current dg is for scalars only
      if (ifflow)  then
         ifield = 1
         call opdssum(vx,vy,vz)
         call opcolv (vx,vy,vz,vmult)
         if (ifsplit) call dsavg(pr)  ! continuous pressure
         if (ifvcor)  call ortho(pr)  ! remove any mean
      endif
 
c     if (ifmhd.and..not.ifdg) then   ! Current dg is for scalars only
      if (ifmhd) then
         ifield = ifldmhd
         call opdssum(bx,by,bz)
         call opcolv (bx,by,bz,vmult)
      endif
 
      if (ifheat.and..not.ifdg) then  ! Don't project if using DG
         ifield = 2
         call dssum(t ,lx1,ly1,lz1)
         call col2 (t ,tmult,ntott)
         do ifield=3,nfield
            if(gsh_fld(ifield).ge.0) then
              call dssum(t(1,1,1,1,ifield-1),lx1,ly1,lz1)
              if(iftmsh(ifield)) then
                call col2 (t(1,1,1,1,ifield-1),tmult,ntott)
              else
                call col2 (t(1,1,1,1,ifield-1),vmult,ntotv)
              endif
            endif
         enddo
      endif
 
c     if (ifpert.and..not.ifdg) then ! Still not DG
      if (ifpert) then
         do jp=1,npert
            ifield = 1
            call opdssum(vxp(1,jp),vyp(1,jp),vzp(1,jp))
            call opcolv (vxp(1,jp),vyp(1,jp),vzp(1,jp),vmult)
 
            if (.not.ifdg) then
               do ifield=2,nfield
                  call dssum(tp(1,ifield-1,jp),lx1,ly1,lz1)
                  if(iftmsh(ifield)) then
                     call col2 (tp(1,ifield-1,jp),tmult,ntott)
                  else
                     call col2 (tp(1,ifield-1,jp),vmult,ntotv)
                  endif
               enddo
            endif
         enddo
      endif
      jp = 0
 
      return
      end