chkptms.f

Go to the documentation of this file.

 
!=======================================================================
      subroutine chkpts_register()
      implicit none
 
      include 'SIZE'
      include 'FRAMELP'
      include 'CHKPTD'
      include 'CHKPTMSD'
 
      ! local variables
      integer lpmid
!-----------------------------------------------------------------------
      ! check if the current module was already registered
      call mntr_mod_is_name_reg(lpmid,chkptms_name)
      if (lpmid.gt.0) then
         call mntr_warn(lpmid,
     $        'module ['//trim(chkptms_name)//'] already registered')
         return
      endif
 
      ! find parent module
      call mntr_mod_is_name_reg(lpmid,chpt_name)
      if (lpmid.le.0) then
        lpmid = 1
        call mntr_abort(lpmid,
     $   'Parent ['//trim(chpt_name)//'] module not registered')
      endif
 
      ! register module
      call mntr_mod_reg(chkptms_id,lpmid,chkptms_name,
     $         'Multi-file checkpointing')
 
      ! register timers
      call mntr_tmr_is_name_reg(lpmid,'CHP_TOT')
      call mntr_tmr_reg(chkptms_tread_id,lpmid,chkptms_id,
     $      'CHP_READ','Checkpointing reading time',.true.)
 
      call mntr_tmr_reg(chkptms_twrite_id,lpmid,chkptms_id,
     $      'CHP_WRITE','Checkpointing writing time',.true.)
 
      ! adjust step delay
      call mntr_set_step_delay(chkptms_snmax)
 
      return
      end subroutine
!=======================================================================
      subroutine chkpts_init
      implicit none
 
      include 'SIZE'            ! NID, NPERT
      include 'TSTEP'           ! ISTEP, NSTEPS
      include 'INPUT'           ! IFPERT, PARAM
      include 'FRAMELP'
      include 'CHKPTD'
      include 'CHKPTMSD'
!-----------------------------------------------------------------------
      ! check if the module was already initialised
      if (chkptms_ifinit) then
         call mntr_warn(chkptms_id,
     $        'module ['//trim(chkptms_name)//'] already initiaised.')
         return
      endif
 
      ! get number of snapshots in a set
      if (param(27).lt.0) then
         chkptms_nsnap = nbdinp
      else
         chkptms_nsnap = chkptms_snmax
      endif
 
      ! we support only one perturbation
      if (ifpert) then
         if (npert.gt.1) call mntr_abort(chkptms_id,
     $         'only single perturbation supported')
      endif
 
      chkptms_ifinit = .true.
 
      return
      end subroutine
!=======================================================================
      logical function chkpts_is_initialised()
      implicit none
 
      include 'SIZE'
      include 'CHKPTMSD'
!-----------------------------------------------------------------------
      chkpts_is_initialised = chkptms_ifinit
 
      return
      end function
!=======================================================================
      subroutine chkpts_write()
      implicit none
 
      include 'SIZE'            !
      include 'TSTEP'           ! ISTEP, NSTEPS
      include 'INPUT'           ! IFMVBD, IFREGUO
      include 'FRAMELP'
      include 'CHKPTD'
      include 'CHKPTMSD'
 
      ! local variables
      integer il, ifile, fnum
      real ltim
      character*132 fname(chkptms_fmax)
      logical ifcoord
      logical ifreguol
 
      character*2 str
      character*200 lstring
 
      integer icalldl
      save    icalldl
      data    icalldl  /0/
 
      ! functions
      real dnekclock
!-----------------------------------------------------------------------
      ! no regular mesh
      ifreguol= ifreguo
      ifreguo = .false.
 
      ! do we write a snapshot
      ifile = 0
      if (chpt_stepc.gt.0.and.chpt_stepc.le.chkptms_nsnap) then
      ! timing
         ltim = dnekclock()
 
         ! file number
         ifile = chkptms_nsnap - chpt_stepc +1
         if (ifile.eq.1) call mntr_log(chkptms_id,lp_inf,
     $                             'Writing checkpoint snapshot')
 
         ! initialise I/O data
         call io_init
 
         ! get set of file names in the snapshot
         call chkptms_set_name(fname, fnum, chpt_set_o, ifile)
 
         ! do we wtrite coordinates; we save coordinates in DNS files only
         if (ifmvbd) then  ! moving boundaries save in every file
            ifcoord = .true.
         elseif (ifile.eq.1) then
            ! perturbation mode with constant base flow - only 1 rsX written
            if (ifpert.and.(.not.ifbase)) then
               if (icalldl.eq.0.and.(.not.chpt_ifrst)) then
                  icalldl = 1
                  ifcoord = .true.
               else
                  call chcopy (fname(1),fname(fnum),132)
                  fnum = 1
                  ifcoord = .false.
               endif
            ! DNS, MHD, perturbation with changing base flow - every first rsX file in snapshot
            else
               ifcoord = .true.
            endif
         else
            if (ifpert.and.(.not.ifbase)) then
               call chcopy (fname(1),fname(fnum),132)
               fnum = 1
            endif
            ifcoord = .false.
         endif
 
         ! write down files
         call chkptms_restart_write(fname, fnum, ifcoord)
 
         ! update output set number
         ! we do it after the last file in the set was sucsesfully written
         if (ifile.eq.chkptms_nsnap) then
            write(str,'(I2)') chpt_set_o+1
            lstring = 'Written checkpoint snapshot number: '//trim(str)
            call mntr_log(chkptms_id,lp_prd,lstring)
         endif
 
         ! timing
         ltim = dnekclock() - ltim
         call mntr_tmr_add(chkptms_twrite_id,1,ltim)
      endif
 
      ! put parameters back
      ifreguo = ifreguol
 
      return
      end subroutine
!=======================================================================
      subroutine chkpts_read()
      implicit none
 
      include 'SIZE'            !
      include 'TSTEP'           ! ISTEP, IF_FULL_PRES
      include 'INPUT'           ! IFREGUO, INITC
      include 'FRAMELP'
      include 'CHKPTD'
      include 'CHKPTMSD'
 
      ! local variables
      integer ifile, fnum, fnuml, il
      real dtratio, epsl, ltim
      parameter(epsl = 0.0001)
      logical ifreguol
      character*132 fname(chkptms_fmax),fnamel(chkptms_fmax)
      character*200 lstring
      integer icalld
      save icalld
      data icalld /0/
 
      !functions
      real dnekclock
!-----------------------------------------------------------------------
      ! no regular mesh; important for file name generation
      ifreguol= ifreguo
      ifreguo = .false.
 
      ! this is multi step restart so check for timestep consistency is necessary
      ! this routine gets the information of pressure mesh as well
      if (chpt_ifrst.and.icalld.eq.0) then
         call chkptms_dt_get
         icalld = 1
      endif
 
      if (chpt_ifrst.and.(istep.lt.chkptms_nsnap)) then
 
         ! timing
         ltim = dnekclock()
 
         ifile = istep+1  ! snapshot number
 
         ! initialise I/O data
         call io_init
 
         ! get set of file names in the snapshot
         call chkptms_set_name(fname, fnum, chpt_set_i, ifile)
 
         ! perturbation mode with constant base flow - only 1 rsX written
         if (ifpert.and.(.not.ifbase)) then
            if (ifile.eq.1) then
               il = 0
               call chkptms_set_name(fnamel, fnuml, il, ifile)
               call chcopy (fname(1),fnamel(1),132)
               fnum = 2
            else
               call chcopy (fname(1),fname(fnum),132)
               fnum = 1
            endif
         endif
 
         call chkptms_restart_read(fname, fnum)
 
         ! check time step consistency
         if(ifile.gt.1.and.chkptms_dtstep(ifile).gt.0.0) then
            dtratio = abs(dt-chkptms_dtstep(ifile))
     $              /chkptms_dtstep(ifile)
            if (dtratio.gt.epsl) then
                write(lstring,*) 'Time step inconsistent, new=',
     $            dt,', old=',chkptms_dtstep(ifile)
               call mntr_warn(chkptms_id,lstring)
              ! possible place to exit if this should be trerated as error
            endif
         endif
 
         ! timing
         ltim = dnekclock() - ltim
         call mntr_tmr_add(chkptms_tread_id,1,ltim)
      endif
 
      ! put parameters back
      ifreguo = ifreguol
      if_full_pres=.false.
 
      return
      end subroutine
!=======================================================================
      subroutine chkptms_dt_get
      implicit none
 
      include 'SIZE'
      include 'INPUT'
      include 'RESTART'
      include 'TSTEP'
      include 'FRAMELP'
      include 'CHKPTD'
      include 'CHKPTMSD'
 
      ! local variables
      integer ifile, ierr
      real timerl(chkptms_snmax), p0thr
      character*132 fname, header
      character*3 prefix
      character*4 dummy
!-----------------------------------------------------------------------
      ! which set of files should be used
      if (ifpert.and.(.not.ifbase)) then
         prefix(1:2)='rp'
      else
         prefix(1:2)='rs'
      endif
 
      ! initialise I/O data
      call io_init
 
      ! collect simulation time from file headers
      do ifile=1,chkptms_nsnap
         call chkpt_fname(fname, prefix, chpt_set_i, ifile, ierr)
         call mntr_check_abort(chkptms_id,ierr,
     $    'dt get; file name error')
 
         ierr = 0
         if (nid.eq.pid00) then
            ! open file
            call addfid(fname,fid0)
            ! add ending character; required by C
            fname = trim(fname)//char(0)
            call byte_open(fname,ierr)
            ! read header
            if (ierr.eq.0) then
               call blank     (header,iheadersize)
               call byte_read (header,iheadersize/4,ierr)
            endif
            ! close the file
            if (ierr.eq.0) call byte_close(ierr)
         endif
         call mntr_check_abort(chkptms_id,ierr,
     $       'dt get; error reading header')
 
         call bcast(header,iheadersize)
         ierr = 0
         if (index(header,'#std').eq.1) then
            read(header,*,iostat=ierr) dummy
     $         ,  wdsizr,nxr,nyr,nzr,nelr,nelgr,timer,istpr
     $         ,  ifiler,nfiler
     $         ,  rdcode      ! 74+20=94
     $         ,  p0thr, chkptms_if_pmesh
         else
            ierr = 1
         endif
         call mntr_check_abort(chkptms_id,ierr,
     $       'dt get; error extracting timer')
 
         timerl(ifile) = timer
      enddo
 
      ! get dt
      do ifile=2,chkptms_nsnap
         chkptms_dtstep(ifile) = timerl(ifile) - timerl(ifile-1)
      enddo
 
      return
      end subroutine
!=======================================================================
      subroutine chkptms_set_name(fname, fnum, nset, ifile)
      implicit none
 
      include 'SIZE'            ! NIO
      include 'INPUT'           ! IFMHD, IFPERT, IFBASE
      include 'FRAMELP'
      include 'CHKPTMSD'
 
      ! argument list
      character*132 fname(chkptms_fmax)
      integer fnum, nset, ifile
 
      ! local variables
      integer ifilel, ierr
 
      character*132 bname
      character*3 prefix
!-----------------------------------------------------------------------
      ! fill fname array with 'rsX' (DNS), 'rpX' (pert.) and 'rbX' (MHD) file names
      if (ifmhd) then
         ! file number
         fnum = 2
 
         ! prefix and name for fluid (DNS)
         prefix(1:2)='rs'
         call chkpt_fname(fname(1), prefix, nset, ifile, ierr)
         call mntr_check_abort(chkptms_id,ierr,
     $        'chkptms_set_name; DNS file name error')
 
         ! prefix and name for magnetic field (MHD)
         prefix(1:2)='rb'
         call chkpt_fname(fname(2), prefix, nset, ifile, ierr)
         call mntr_check_abort(chkptms_id,ierr,
     $        'chkptms_set_name; MHD file name error')
 
      elseif (ifpert) then
         ! file number
         ! I assume only single perturbation
         fnum = 2
 
         ! prefix and name for base flow (DNS)
         prefix(1:2)='rs'
         if (ifbase) then
            ifilel = ifile
         else
            ifilel =1
         endif
         call chkpt_fname(fname(1), prefix, nset, ifilel, ierr)
         call mntr_check_abort(chkptms_id,ierr,
     $        'chkptms_set_name; base flow file name error')
 
         ! prefix and name for perturbation
         prefix(1:2)='rp'
         call chkpt_fname(fname(2), prefix, nset, ifile, ierr)
         call mntr_check_abort(chkptms_id,ierr,
     $        'chkptms_set_name; perturbation file name error')
 
      else                ! DNS
         fnum = 1
 
         ! create prefix and name for DNS
         prefix(1:2)='rs'
         call chkpt_fname(fname(1), prefix, nset, ifile, ierr)
         call mntr_check_abort(chkptms_id,ierr,
     $        'chkptms_set_name; DNS file name error')
 
      endif
 
      return
      end subroutine
!=======================================================================
      subroutine chkpt_fname(fname, prefix, nset, ifile, ierr)
      implicit none
 
      include 'SIZE'            ! NIO
      include 'INPUT'           ! SESSION
      include 'CHKPTD'
      include 'CHKPTMSD'
 
      ! argument list
      character*132 fname
      character*3 prefix
      integer nset, ifile, ierr
 
      ! local variables
      character*132 bname    ! base name
      character*132 fnamel   ! local file name
      character*3 prefixl    ! local prefix
      integer itmp
 
      character*6  str
 
      character*(*) kst
      parameter(kst='0123456789abcdefx')
!-----------------------------------------------------------------------
      ! create prefix and name for DNS
      ierr = 0
      prefixl(1:2) = prefix(1:2)
      itmp=min(17,chpt_nset*chkptms_nsnap) + 1
      prefixl(3:3)=kst(itmp:itmp)
 
      ! get base name (SESSION)
      bname = trim(adjustl(session))
 
      call io_mfo_fname(fnamel,bname,prefixl,ierr)
      if (ierr.ne.0) then
         call mntr_error(chkptms_id,'chkpt_fname; file name error')
         return
      endif
 
      write(str,'(i5.5)') chkptms_nsnap*nset+ifile
      fname = trim(fnamel)//trim(str)
 
      return
      end subroutine
!=======================================================================
      subroutine chkptms_restart_write(fname, fnum, ifcoord)
      implicit none
 
      include 'SIZE'
      include 'RESTART'
      include 'TSTEP'
      include 'INPUT'
      include 'CHKPTMSD'
 
      ! argument list
      character*132 fname(chkptms_fmax)
      integer fnum
      logical ifcoord
 
      ! local variables
      integer lwdsizo
      integer ipert, il
      ! which set of variables do we write: DNS (1), MHD (2) or perturbation (3)
      integer chktype
 
      logical lif_full_pres, lifxyo, lifpo, lifvo, lifto,
     $        lifpsco(LDIMT1)
!-----------------------------------------------------------------------
      ! adjust I/O parameters
      lwdsizo = wdsizo
      wdsizo  = 8
      lif_full_pres = if_full_pres
      if_full_pres = .true.
      lifxyo = ifxyo
      lifpo= ifpo
      ifpo = .true.
      lifvo= ifvo
      ifvo = .true.
      lifto= ifto
      ifto = ifheat
      do il=1,npscal
         lifpsco(il)= ifpsco(il)
         ifpsco(il) = .true.
      enddo
 
      if (ifmhd) then
         ! DNS first
         ifxyo = ifcoord
         chktype = 1
         call chkptms_mfo(fname(1),chktype,ipert)
 
         ! MHD
         ifxyo = .false.
         chktype = 2
         call chkptms_mfo(fname(2),chktype,ipert)
 
      elseif (ifpert) then
         ! DNS first
         if (fnum.eq.2) then
            ifxyo = ifcoord
            chktype = 1
            call chkptms_mfo(fname(1),chktype,ipert)
         endif
 
         ! perturbation
         ifxyo = .false.
         chktype = 3
         ipert = 1
         call chkptms_mfo(fname(fnum),chktype,ipert)
      else ! DNS
         ! write only one set of files
         if (fnum.ne.1) call mntr_abort(chkptms_id,
     $        'chkptms_restart_write; too meny files for DNS')
 
         ifxyo = ifcoord
         chktype = 1
         call chkptms_mfo(fname(1),chktype,ipert)
      endif
 
      ! restore I/O parameters
      wdsizo = lwdsizo
      if_full_pres = lif_full_pres
      ifxyo = lifxyo
      ifpo = lifpo
      ifvo = lifvo
      ifto = lifto
      do il=1,npscal
         ifpsco(il) = lifpsco(il)
      enddo
 
      return
      end subroutine
!=======================================================================
      subroutine chkptms_restart_read(fname, fnum)
      implicit none
 
      include 'SIZE'
      include 'RESTART'
      include 'TSTEP'
      include 'INPUT'
      include 'CHKPTMSD'
 
      ! argument list
      character*132 fname(chkptms_fmax)
      integer fnum
 
      ! local variables
      integer ndumps, ipert, il
      ! which set of variables do we write: DNS (1), MHD (2) or perturbation (3)
      integer chktype
      character*132 fnamel
!-----------------------------------------------------------------------
      if (ifmhd) then
         ! DNS first
         chktype = 1
         call sioflag(ndumps,fnamel,fname(1))
         call chkptms_mfi(fnamel,chktype,ipert)
 
         ! MHD
         chktype = 2
         call sioflag(ndumps,fnamel,fname(2))
         call chkptms_mfi(fnamel,chktype,ipert)
 
      elseif (ifpert) then
         ! DNS first
         if (fnum.eq.2) then
            chktype = 1
            call sioflag(ndumps,fnamel,fname(1))
            call chkptms_mfi(fnamel,chktype,ipert)
         endif
 
         ! perturbation
         chktype = 3
         ipert = 1
         call sioflag(ndumps,fnamel,fname(fnum))
         call chkptms_mfi(fnamel,chktype,ipert)
      else ! DNS
         ! read only one set of files
         if (fnum.ne.1) call mntr_abort(chkptms_id,
     $        'chkptms_restart_read; too meny files for DNS')
 
         chktype = 1
         call sioflag(ndumps,fnamel,fname(1))
         call chkptms_mfi(fnamel,chktype,ipert)
      endif
 
      return
      end subroutine
!=======================================================================
      subroutine chkptms_mfo(fname,chktype,ipert)
      implicit none
 
      include 'SIZE'
      include 'INPUT'
      include 'PARALLEL'
      include 'RESTART'
      include 'TSTEP'
      include 'GEOM'
      include 'SOLN'
      include 'FRAMELP'
      include 'CHKPTMSD'
 
      ! argumnt list
      character*132 fname
      integer chktype, ipert
 
      ! local variables
      integer ierr, il, itmp
      integer ioflds
      integer*8 offs,nbyte
      real dnbyte, tiostart, tio
 
      ! functions
      real dnekclock_sync, glsum
 
      real pm1(lx1,ly1,lz1,lelv)
      common /scrcg/ pm1
!-----------------------------------------------------------------------
      ! simple timing
      tiostart=dnekclock_sync()
 
      ! set elelemnt size
      nxo  = nx1
      nyo  = ny1
      nzo  = nz1
 
      ! open file
      call io_mbyte_open(fname,ierr)
      call mntr_check_abort(chkptms_id,ierr,
     $   'chkptms_mfo; file not opened')
 
      ! write a header and create element mapping
      call mfo_write_hdr
 
      ! set header offset
      offs = iheadersize + 4 + isize*nelgt
      ioflds = 0
 
      ! write fields
      ! coordinates
      if (ifxyo) then
         call io_mfov(offs,xm1,ym1,zm1,nx1,ny1,nz1,nelt,nelgt,ndim)
         ioflds = ioflds + ndim
      endif
 
      ! velocity, magnetic field, perturbation velocity
      if (ifvo) then
         if (chktype.eq.1) then
            call io_mfov(offs,vx,vy,vz,nx1,ny1,nz1,nelt,nelgt,ndim)
         elseif(chktype.eq.2) then
            call io_mfov(offs,bx,by,bz,nx1,ny1,nz1,nelt,nelgt,ndim)
         elseif(chktype.eq.3) then
            call io_mfov(offs,vxp(1,ipert),vyp(1,ipert),
     $                       vzp(1,ipert),nx1,ny1,nz1,nelt,nelgt,ndim)
         endif
         ioflds = ioflds + ndim
      endif
 
      ! pressure
      if (ifpo) then
         if (chktype.eq.1) then
            ! copy array if necessary
            if (ifsplit) then
               itmp = nx1*ny1*nz1*lelv
               call io_mfos(offs,pr,nx2,ny2,nz2,nelt,nelgt,ndim)
            else
               itmp = nx1*ny1*nz1*lelv
               call rzero(pm1,itmp)
               itmp = nx2*ny2*nz2
               do il=1,nelv
                  call copy(pm1(1,1,1,il),pr(1,1,1,il),itmp)
               enddo
               call io_mfos(offs,pm1,nx1,ny1,nz1,nelt,nelgt,ndim)
            endif
         elseif(chktype.eq.2) then
            ! copy array
            itmp = nx1*ny1*nz1*lelv
            call rzero(pm1,itmp)
            itmp = nx2*ny2*nz2
            do il=1,nelv
               call copy(pm1(1,1,1,il),pm(1,1,1,il),itmp)
            enddo
            call io_mfos(offs,pm1,nx1,ny1,nz1,nelt,nelgt,ndim)
         elseif(chktype.eq.3) then
            ! copy array
            itmp = nx1*ny1*nz1*lelv
            call rzero(pm1,itmp)
            itmp = nx2*ny2*nz2
            do il=1,nelv
               call copy(pm1(1,1,1,il),prp(1+itmp*(il-1),ipert),itmp)
            enddo
            call io_mfos(offs,pm1,nx1,ny1,nz1,nelt,nelgt,ndim)
         endif
         ioflds = ioflds + 1
      endif
 
      if (chktype.ne.2) then
         if (ifto) then
            if (chktype.eq.1) then
               call io_mfos(offs,t,nx1,ny1,nz1,nelt,nelgt,ndim)
            elseif(chktype.eq.3) then
               call io_mfos(offs,tp(1,1,ipert),
     $                          nx1,ny1,nz1,nelt,nelgt,ndim)
            endif
            ioflds = ioflds + 1
         endif
 
         do il=1,ldimt-1
            if (ifpsco(il)) then
               if (chktype.eq.1) then
                  call io_mfos(offs,t(1,1,1,1,il+1),
     $                 nx1,ny1,nz1,nelt,nelgt,ndim)
               elseif(chktype.eq.3) then
                  call io_mfos(offs,tp(1,il+1,ipert),
     $                 nx1,ny1,nz1,nelt,nelgt,ndim)
               endif
               ioflds = ioflds + 1
            endif
         enddo
      endif
      dnbyte = 1.*ioflds*nelt*wdsizo*nx1*ny1*nz1
 
      ! possible place for metadata
 
      ! close file
      call io_mbyte_close(ierr)
      call mntr_check_abort(chkptms_id,ierr,
     $    'chkptms_mfo; file not closed')
 
      ! stamp the log
      tio = dnekclock_sync()-tiostart
      if (tio.le.0) tio=1.
 
      dnbyte = glsum(dnbyte,1)
      dnbyte = dnbyte + iheadersize + 4. + isize*nelgt
      dnbyte = dnbyte/1024/1024
 
      call mntr_log(chkptms_id,lp_prd,'Checkpoint written:')
      call mntr_logr(chkptms_id,lp_vrb,'file size (MB) = ',dnbyte)
      call mntr_logr(chkptms_id,lp_vrb,'avg data-throughput (MB/s) = ',
     $     dnbyte/tio)
      call mntr_logi(chkptms_id,lp_vrb,'io-nodes = ',nfileo)
 
      return
      end subroutine
!=======================================================================
      subroutine chkptms_mfi(fname,chktype,ipert)
      implicit none
 
      include 'SIZE'
      include 'INPUT'
      include 'PARALLEL'
      include 'RESTART'
      include 'TSTEP'
      include 'GEOM'
      include 'SOLN'
      include 'FRAMELP'
      include 'CHKPTMSD'
 
      ! argumnt list
      character*132 fname
      integer chktype, ipert
 
      ! local variables
      integer ierr, il, jl, kl
      integer itmp1, itmp2, itmp3
      integer ioflds
      integer*8 offs0,offs,nbyte
      real dnbyte, tiostart, tio
      logical ifskip
 
      !     scratch arrays
      integer lwkv
      parameter(lwkv = lx1*ly1*lz1*lelt)
      real wkv1(lwkv),wkv2(lwkv),wkv3(lwkv)
      common /scruz/ wkv1,wkv2,wkv3
 
      ! functions
      real dnekclock_sync, glsum
!-----------------------------------------------------------------------
      ! simple timing
      tiostart=dnekclock_sync()
 
      ! open file and read header; some operations related to header are
      ! performed in chkpts_read; it is not the optimal way, but I don't want
      ! to modify mfi_prepare
      call mfi_prepare(fname)
 
      ! set header offset
      offs = iheadersize + 4 + isize*nelgr
      ioflds = 0
 
      ! read fields
      ! coordinates
      if (ifgetxr) then
         ifskip = .not.ifgetx
         ! skip coordinates if you need to interpolate them;
         ! I assume current coordinates are more exact than the interpolated one
         if ((nxr.eq.lx1).and.(nyr.eq.ly1).and.(nzr.eq.lz1)) then
            call io_mfiv(offs,xm1,ym1,zm1,lx1,ly1,lz1,lelt,ifskip)
         else
            ifskip=.true.
            call io_mfiv(offs,wkv1,wkv2,wkv3,nxr,nyr,nzr,lelt,ifskip)
         endif
         ioflds = ioflds + ldim
      endif
 
      ! velocity, magnetic field, perturbation velocity
      if (ifgetur) then
         ifskip = .not.ifgetu
         if ((nxr.eq.lx1).and.(nyr.eq.ly1).and.(nzr.eq.lz1)) then
            ! unchanged resolution
            ! read field directly to the variables
            if (chktype.eq.1) then
               call io_mfiv(offs,vx,vy,vz,lx1,ly1,lz1,lelv,ifskip)
            elseif(chktype.eq.2) then
               call io_mfiv(offs,bx,by,bz,lx1,ly1,lz1,lelv,ifskip)
            elseif(chktype.eq.3) then
               call io_mfiv(offs,vxp(1,ipert),vyp(1,ipert),vzp(1,ipert)
     $                   ,lx1,ly1,lz1,lelv,ifskip)
            endif
         else
            ! modified resolution
            ! read field to tmp array
            call io_mfiv(offs,wkv1,wkv2,wkv3,nxr,nyr,nzr,lelt,ifskip)
 
            ! interpolate
            if (ifgetu) then
               if (chktype.eq.1) then
                  call chkptms_map_gll(vx,wkv1,nxr,nzr,nelv)
                  call chkptms_map_gll(vy,wkv2,nxr,nzr,nelv)
                  if (if3d) call chkptms_map_gll(vz,wkv3,nxr,nzr,nelv)
               elseif(chktype.eq.2) then
                  call chkptms_map_gll(bx,wkv1,nxr,nzr,nelv)
                  call chkptms_map_gll(by,wkv2,nxr,nzr,nelv)
                  if (if3d) call chkptms_map_gll(bz,wkv3,nxr,nzr,nelv)
               elseif(chktype.eq.3) then
                  call chkptms_map_gll(vxp(1,ipert),wkv1,nxr,nzr,nelv)
                  call chkptms_map_gll(vyp(1,ipert),wkv2,nxr,nzr,nelv)
                  if (if3d) call chkptms_map_gll(vzp(1,ipert),wkv3,
     $                        nxr,nzr,nelv)
               endif
            endif
         endif
         ioflds = ioflds + ndim
      endif
 
      ! pressure
      if (ifgetpr) then
         ifskip = .not.ifgetp
 
         if (chkptms_if_pmesh) then
            ! file contains pressure on nxr-2 (GL) mesh
            ! read field to tmp array
            call io_mfis(offs,wkv1,nxr,nyr,nzr,lelt,ifskip)
 
            if (ifgetp) then
               if ((nxr.eq.lx1).and.(nyr.eq.ly1).and.(nzr.eq.lz1)) then
                  ! unchanged resolution
                  if (ifsplit) then
                     if (chktype.eq.1) then
                        !interpolate GL to GLL
                        itmp1 = nx1*ny1*nz1
                        jl = 1
                        do il=1,nelv
                           call map21t (pr(1,1,1,il),wkv1(jl),il)
                           jl = jl +itmp1
                        enddo
                     endif
                  else
                     ! remove zeros
                     itmp1 = nx1*ny1*nz1
                     itmp2 = nx2*ny2*nz2
                     jl = 1
                     if(chktype.eq.1) then
                        do il=1,nelv
                           call copy(pr(1,1,1,il),wkv1(jl),itmp2)
                           jl = jl + itmp1
                        enddo
                     elseif(chktype.eq.2) then
                        do il=1,nelv
                           call copy(pm(1,1,1,il),wkv1(jl),itmp2)
                           jl = jl + itmp1
                        enddo
                     elseif(chktype.eq.3) then
                        kl = 1
                        do il=1,nelv
                           call copy(prp(kl,ipert),wkv1(jl),itmp2)
                           jl = jl + itmp1
                           kl = kl + itmp2
                        enddo
                     endif
                  endif
               else
                  ! modified resolution
                  ! remove zeros
                  itmp1 = nxr*nyr*nzr
                  itmp3 = max(nzr-2,1)
                  itmp2 = (nxr-2)*(nyr-2)*itmp3
                  jl = 1
                  kl = 1
                  do il=1,nelv
                     call copy(wkv2(kl),wkv1(jl),itmp2)
                     jl = jl + itmp1
                     kl = kl + itmp2
                  enddo
 
                  if (ifsplit) then
                     if (chktype.eq.1) then
                        ! interpolate on GL mesh
                        call chkptms_map_gl(wkv1,wkv2,nxr-2,itmp3,nelv)
 
                        !interpolate GL to GLL
                        itmp2 = nx2*ny2*nz2
                        jl = 1
                        do il=1,nelv
                           call map21t (pr(1,1,1,il),wkv1(jl),il)
                           jl = jl +itmp2
                        enddo
                     endif
                  else
                     ! interpolate on GL mesh
                     if (chktype.eq.1) then
                        call chkptms_map_gl(pr,wkv2,nxr-2,itmp3,nelv)
                     elseif(chktype.eq.2) then
                        call chkptms_map_gl(pm,wkv2,nxr-2,itmp3,nelv)
                     elseif(chktype.eq.3) then
                        call chkptms_map_gl(prp(1,ipert),wkv2,nxr-2,
     $                       itmp3,nelv)
                     endif
                  endif
               endif
            endif
 
         else
            ! file contains pressure on nxr (GLL) mesh
            if ((nxr.eq.lx1).and.(nyr.eq.ly1).and.(nzr.eq.lz1)) then
               ! unchanged resolution
               if (ifsplit) then
                  ! read field directly to the variables
                  if (chktype.eq.1) then
                     call io_mfis(offs,pr,lx1,ly1,lz1,lelv,ifskip)
                  endif
               else
                  ! read field to tmp array
                  call io_mfis(offs,wkv1,nxr,nyr,nzr,lelt,ifskip)
 
                  ! interpolate GLL to GL
                  if (ifgetp) then
                     itmp2 = nx1*ny1*nz1
                     itmp2 = nx2*ny2*nz2
                     jl = 1
                     if (chktype.eq.1) then
                        do il=1,nelv
                           call map12 (pr(1,1,1,il),wkv1(jl),il)
                           jl = jl +itmp1
                        enddo
                     elseif(chktype.eq.2) then
                        do il=1,nelv
                           call map12 (pm(1,1,1,il),wkv1(jl),il)
                           jl = jl +itmp1
                        enddo
                     elseif(chktype.eq.3) then
                        kl = 1
                        do il=1,nelv
                           call map12 (prp(kl,ipert),wkv1(jl),il)
                           jl = jl +itmp1
                           kl = kl +itmp2
                        enddo
                     endif
                  endif
               endif
 
            else
               ! modified resolution
               ! read field to tmp array
               call io_mfis(offs,wkv1,nxr,nyr,nzr,lelt,ifskip)
 
               if (ifgetp) then
                  if (ifsplit) then
                     ! interpolate on GLL
                     call chkptms_map_gll(pr,wkv1,nxr,nzr,nelv)
                  else
                     ! interpolate on GLL
                     call chkptms_map_gll(wkv2,wkv1,nxr,nzr,nelv)
 
                     ! interpolate GLL to GL
                     itmp2 = nx1*ny1*nz1
                     itmp2 = nx2*ny2*nz2
                     jl = 1
                     if (chktype.eq.1) then
                        do il=1,nelv
                           call map12 (pr(1,1,1,il),wkv2(jl),il)
                           jl = jl +itmp1
                        enddo
                     elseif(chktype.eq.2) then
                        do il=1,nelv
                           call map12 (pm(1,1,1,il),wkv2(jl),il)
                           jl = jl +itmp1
                        enddo
                     elseif(chktype.eq.3) then
                        kl = 1
                        do il=1,nelv
                           call map12 (prp(kl,ipert),wkv2(jl),il)
                           jl = jl +itmp1
                           kl = kl +itmp2
                        enddo
                     endif
                  endif
               endif
            endif
         endif
 
         ioflds = ioflds + 1
      endif
 
      if (chktype.ne.2) then
         if (ifgettr) then
            ifskip = .not.ifgett
            if ((nxr.eq.lx1).and.(nyr.eq.ly1).and.(nzr.eq.lz1)) then
               ! unchanged resolution
               ! read field directly to the variables
               if (chktype.eq.1) then
                  call io_mfis(offs,t,lx1,ly1,lz1,lelt,ifskip)
               elseif(chktype.eq.3) then
                  call io_mfis(offs,tp(1,1,ipert),lx1,ly1,lz1,lelt,
     $                 ifskip)
               endif
            else
               ! modified resolution
               ! read field to tmp array
               call io_mfis(offs,wkv1,nxr,nyr,nzr,lelt,ifskip)
 
               ! interpolate
               if (ifgett) then
                  if (chktype.eq.1) then
                     call chkptms_map_gll(t,wkv1,nxr,nzr,nelt)
                  elseif(chktype.eq.3) then
                     call chkptms_map_gll(tp(1,1,ipert),wkv1,nxr,nzr,
     $               nelt)
                  endif
               endif
            endif
            ioflds = ioflds + 1
         endif
 
         do il=1,ldimt-1
            if (ifgtpsr(il)) then
               ifskip = .not.ifgtps(il)
               if ((nxr.eq.lx1).and.(nyr.eq.ly1).and.(nzr.eq.lz1)) then
                  ! unchanged resolution
                  ! read field directly to the variables
                  if (chktype.eq.1) then
                     call io_mfis(offs,t(1,1,1,1,il+1),lx1,ly1,lz1,
     $                            lelt,ifskip)
                  elseif(chktype.eq.3) then
                     call io_mfis(offs,tp(1,il+1,ipert),lx1,ly1,lz1,
     $                            lelt,ifskip)
                  endif
               else
                  ! modified resolution
                  ! read field to tmp array
                  call io_mfis(offs,wkv1,nxr,nyr,nzr,lelt,ifskip)
 
                  ! interpolate
                  if (ifgtps(il)) then
                     if (chktype.eq.1) then
                        call chkptms_map_gll(t(1,1,1,1,il+1),wkv1,
     $                                    nxr,nzr,nelt)
                     elseif(chktype.eq.3) then
                        call chkptms_map_gll(tp(1,il+1,ipert),wkv1,
     $                                    nxr,nzr,nelt)
                     endif
                  endif
               endif
               ioflds = ioflds + 1
            endif
         enddo
      endif
 
      if (ifgtim) time = timer
 
      ! close file
      call io_mbyte_close(ierr)
      call mntr_check_abort(chkptms_id,ierr,
     $     'chkptms_mfi; file not closed')
 
      ! stamp the log
      tio = dnekclock_sync()-tiostart
      if (tio.le.0) tio=1.
 
      if(nid.eq.pid0r) then
         dnbyte = 1.*ioflds*nelr*wdsizr*nxr*nyr*nzr
      else
         dnbyte = 0.0
      endif
 
      dnbyte = glsum(dnbyte,1)
      dnbyte = dnbyte + iheadersize + 4. + isize*nelgt
      dnbyte = dnbyte/1024/1024
 
      call mntr_log(chkptms_id,lp_prd,'Checkpoint read:')
      call mntr_logr(chkptms_id,lp_vrb,'avg data-throughput (MB/s) = ',
     $     dnbyte/tio)
      call mntr_logi(chkptms_id,lp_vrb,'io-nodes = ',nfileo)
 
      if (ifaxis) call chkptms_axis_interp_ic()
 
      return
      end subroutine
!=======================================================================
      subroutine chkptms_map_gll(xf,yf,nxr,nzr,nel)
      implicit none
 
      include 'SIZE'
      include 'INPUT'
      include 'IXYZ'
      include 'WZ'
 
      ! argument lists
      integer nxr, nzr, nel
      real xf(lx1,ly1,lz1,nel), yf(nxr,nxr,nzr,nel)
 
      ! local variables
      integer nyzr, nxy2
      integer ie, iz, izoff
 
      ! moddificaion flag
      integer nold
      save    nold
      data    nold /0/
 
      ! work arrays
      integer lxr, lyr, lzr, lxyzr
      parameter(lxr=lx1+6)
      parameter(lyr=ly1+6)
      parameter(lzr=lz1+6)
      parameter(lxyzr=lxr*lyr*lzr)
      real txa(lxyzr),txb(lx1,ly1,lzr),zgmr(lxr),wgtr(lxr)
      common /ctmp0/  txa, txb, zgmr, wgtr
 
      ! interpolation arrays
      real ires(lxr*lxr)  ,itres(lxr*lxr)
      common /ctmpabm1/ ires, itres
!-----------------------------------------------------------------------
      nyzr = nxr*nzr
      nxy2 = lx1*ly1
 
      if (nxr.ne.nold) then
         nold=nxr
         call zwgll(zgmr,wgtr,nxr)
         call igllm(ires,itres,zgmr,zgm1,nxr,lx1,nxr,lx1)
      endif
 
      do ie=1,nel
         call mxm (ires,lx1,yf(1,1,1,ie),nxr,txa,nyzr)
         do iz=1,nzr
            izoff = 1 + (iz-1)*lx1*nxr
            call mxm (txa(izoff),lx1,itres,nxr,txb(1,1,iz),ly1)
         enddo
         if (if3d) then
            call mxm (txb,nxy2,itres,nzr,xf(1,1,1,ie),lz1)
         else
            call copy(xf(1,1,1,ie),txb,nxy2)
         endif
      enddo
 
      return
      end subroutine
!=======================================================================
      subroutine chkptms_map_gl(xf,yf,nxr,nzr,nel)
      implicit none
 
      include 'SIZE'
      include 'INPUT'
      include 'IXYZ'
      include 'WZ'
 
      ! argument lists
      integer nxr, nzr, nel
      real xf(lx2,ly2,lz2,nel), yf(nxr,nxr,nzr,nel)
 
      ! local variables
      integer nyzr, nxy2
      integer ie, iz, izoff
 
      ! moddificaion flag
      integer nold
      save    nold
      data    nold /0/
 
      ! work arrays
      integer lxr, lyr, lzr, lxyzr
      parameter(lxr=lx2+6)
      parameter(lyr=ly2+6)
      parameter(lzr=lz2+6)
      parameter(lxyzr=lxr*lyr*lzr)
      real txa(lxyzr),txb(lx2,ly2,lzr),zgmr(lxr),wgtr(lxr)
      common /ctmp0/  txa, txb, zgmr, wgtr
 
      ! interpolation arrays
      real ires(lxr,lxr)  ,itres(lxr,lxr)
      common /ctmpabm2/ ires, itres
!-----------------------------------------------------------------------
      nyzr = nxr*nzr
      nxy2 = lx2*ly2
 
      if (nxr.ne.nold) then
         nold=nxr
         call zwgl   (zgmr,wgtr,nxr)
         call iglm   (ires,itres,zgmr,zgm2,nxr,lx2,nxr,lx2)
      endif
 
      do ie=1,nel
         call mxm (ires,lx2,yf(1,1,1,ie),nxr,txa,nyzr)
         do iz=1,nzr
            izoff = 1 + (iz-1)*lx2*nxr
            call mxm (txa(izoff),lx2,itres,nxr,txb(1,1,iz),ly2)
         enddo
         if (if3d) then
            call mxm (txb,nxy2,itres,nzr,xf(1,1,1,ie),lz2)
         else
            call copy(xf(1,1,1,ie),txb,nxy2)
         endif
      enddo
 
      return
      end subroutine
!=======================================================================
      subroutine chkptms_axis_interp_ic()
      implicit none
 
      include 'SIZE'
      include 'INPUT'
      include 'TSTEP'
      include 'RESTART'
      include 'SOLN'
      include 'GEOM'
      include 'IXYZ'
 
      ! scratch space
      real axism1 (lx1,ly1), axism2 (lx2,ly2), ialj2 (ly2,ly2),
     $     iatlj2(ly2,ly2), tmp(ly2,ly2)
      common /ctmp0/ axism1, axism2, ialj2, iatlj2, tmp
 
      ! local variables
      integer el, ips, is1
!-----------------------------------------------------------------------
      if (.not.ifaxis) return
 
      ! get interpolation operators between Gauss-Lobatto Jacobi
      ! and and Gauss Legendre poits (this is missing in genwz)
      !call invmt(iajl2 ,ialj2 ,tmp ,ly2)
      call invmt(iatjl2,iatlj2,tmp,ly2)
 
      do el=1,nelv
         if (ifrzer(el)) then
           if (ifgetx) then
             call mxm(xm1(1,1,1,el),nx1,iatlj1,ny1,axism1,ny1)
             call copy(xm1(1,1,1,el),axism1,nx1*ny1)
             call mxm(ym1(1,1,1,el),nx1,iatlj1,ny1,axism1,ny1)
             call copy(ym1(1,1,1,el),axism1,nx1*ny1)
           endif
           if (ifgetu) then
             call mxm(vx(1,1,1,el),nx1,iatlj1,ny1,axism1,ny1)
             call copy(vx(1,1,1,el),axism1,nx1*ny1)
             call mxm(vy(1,1,1,el),nx1,iatlj1,ny1,axism1,ny1)
             call copy(vy(1,1,1,el),axism1,nx1*ny1)
           endif
           if (ifgetw) then
             call mxm(vz(1,1,1,el),nx1,iatlj1,ny1,axism1,ny1)
             call copy(vz(1,1,1,el),axism1,nx1*ny1)
           endif
           if (ifgetp) then
             if (ifsplit) then
                call mxm(pr(1,1,1,el),nx1,iatlj1,ny1,axism1,ny1)
                call copy(pr(1,1,1,el),axism1,nx1*ny1)
             else
                call mxm(pr(1,1,1,el),nx2,iatlj2,ny2,axism2,ny2)
                call copy(pr(1,1,1,el),axism2,nx2*ny2)
             endif
           endif
           if (ifgett) then
             call mxm(t(1,1,1,el,1),nx1,iatlj1,ny1,axism1,ny1)
             call copy(t(1,1,1,el,1),axism1,nx1*ny1)
           endif
           do ips=1,npscal
            is1 = ips + 1
            if (ifgtps(ips)) then
             call mxm(t(1,1,1,el,is1),nx1,iatlj1,ny1,axism1,ny1)
             call copy(t(1,1,1,el,is1),axism1,nx1*ny1)
            endif
           enddo
         endif
      enddo
 
      return
      end subroutine
!=======================================================================