! This file is part of GTC version 4.4 
! GTC version 4.4 is released under the 3-Clause BSD license:

! Copyright (c) 2002,2010,2016, GTC Team (team leader: Zhihong Lin, zhihongl@uci.edu)
! All rights reserved.

! Redistribution and use in source and binary forms, with or without 
! modification, are permitted provided that the following conditions are met:

! 1. Redistributions of source code must retain the above copyright notice, 
!    this list of conditions and the following disclaimer.

! 2. Redistributions in binary form must reproduce the above copyright notice, 
!    this list of conditions and the following disclaimer in the documentation 
!    and/or other materials provided with the distribution.

! 3. Neither the name of the GTC Team nor the names of its contributors may be 
!    used to endorse or promote products derived from this software without 
!    specific prior written permission.
! ==============================================================================

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!the unformatted bin file
subroutine restart_native(iop)
  use system_env
  use global_parameters
  use particle_array
  use field_array
  use conservative_scheme,only:sapara_na,apara_na,phiind
#ifdef ADIOS2
  use adios2_common
#endif
  !use simple_error_handler
  use mpi
  implicit none

  character(*),intent(in)::iop
  character(len=80) :: restart_fname
  character(len=18) cdum

  if(mype < 10)then
     write(cdum,'("DATA_RESTART.0000",i1)')mype
  elseif(mype < 100)then
     write(cdum,'("DATA_RESTART.000",i2)')mype
  elseif(mype < 1000)then
     write(cdum,'("DATA_RESTART.00",i3)')mype
  elseif(mype < 10000)then
     write(cdum,'("DATA_RESTART.0",i4)')mype
  else
     write(cdum,'("DATA_RESTART.",i5)')mype
  endif

!!XY rolling restart: save two copies of restart data
#ifdef ADIOS2
  if(mod(irest,2)==0)then
     restart_fname="restart_dir1/DATA_RESTART.bp"
  else
     restart_fname="restart_dir2/DATA_RESTART.bp"
  endif
#else
  if(mod(irest,2)==0)then
     restart_fname="restart_dir1/"//trim(cdum)
  else
     restart_fname="restart_dir2/"//trim(cdum)
  endif
  call MPI_BARRIER(MPI_COMM_WORLD,ierror)
#endif

! record particle information for future restart run
  if(iop=="write")then
#ifdef ADIOS2
    call restart_adios2_write(restart_fname)
#else
     open(222,file=trim(restart_fname),status='replace',form='unformatted')
     !$acc update host(zion,zonali,zonalci)
     write(222)mi,me,rdtemi,rddeni,pfluxi,phi,phip00,phi00,zonali,zonalci,meshti,meshte,meshni,meshne,pmarki,markeri,markerit,&
          zion(:,1:mi)

     if(nhybrid>0)then
     !$acc update host(zelectron,zonale,zonalce,phisave,dnesave,pressureepara,pressureeperp)
       write(222)rdteme,rddene,pfluxe,zonale,zonalce,phisave,dnesave,zelectron(:,1:me),&
       pmarke,markere,tfrac,tfracn,&
       densitye,pressureepara,pressureeperp
     endif

     if(fload>0)then
     !$acc update host(zfast,zonalf,zonalcf)
       write(222)mf,rdtemf,rddenf,zonalf,zonalcf,meshtf,meshnf,pfluxf,zfast(:,1:mf),&
       pmarkf,markerf
     endif

     if(feload>0)then
     !$acc update host(zfaste,zonalfe,zonalcfe)
       write(222)mfe,rdtemfe,rddenfe,zonalfe,zonalcfe,meshtfe,meshnfe,pfluxfe,zfaste(:,1:mfe),&
       pmarkfe,markerfe,fetfrac,fetfracn
     endif

     if(magnetic>0)write(222)apara,sapara,sdelapara,fluidne,sfluidne,fluidue,deltapsi,sdeltapsi,aparanl,apara00,fluidne00,gradgaugef
     if(deltabpara==1)write(222)bpara,bparasave
     if(cs_method==1)write(222)sapara_na,apara_na,phiind
     close(222)
#endif
     if(mype==0)write(gtcout,*)'write to ',trim(restart_fname)
     call restart_hist

! read particle information to restart previous run
  else
#ifdef ADIOS2
    call restart_adios2_read(restart_fname)
#else
     if(mype==0)write(gtcout,*)'read in ',trim(restart_fname)
     open(333,file=trim(restart_fname),status='old',form='unformatted')
     read(333)mi,me,rdtemi,rddeni,pfluxi,phi,phip00,phi00,zonali,zonalci,meshti,meshte,meshni,meshne,pmarki,markeri,markerit,&
          zion(:,1:mi)
     if(nhybrid>0)then
          read(333)rdteme,rddene,pfluxe,zonale,zonalce,phisave,dnesave,zelectron(:,1:me),&
          pmarke,markere,tfrac,tfracn,&
          densitye,pressureepara,pressureeperp
          !$acc update device(phisave,dnesave,pressureepara,pressureeperp,densitye)
     endif
     if(fload>0)read(333)mf,rdtemf,rddenf,zonalf,zonalcf,meshtf,meshnf,pfluxf,zfast(:,1:mf),&
          pmarkf,markerf
     if(feload>0)read(333)mfe,rdtemfe,rddenfe,zonalfe,zonalcfe,meshtfe,meshnfe,pfluxfe,zfaste(:,1:mfe),&
          pmarkfe,markerfe,fetfrac,fetfracn
     if(magnetic>0)read(333)apara,sapara,sdelapara,fluidne,sfluidne,fluidue,deltapsi,sdeltapsi,aparanl,apara00,fluidne00,gradgaugef
     if(deltabpara==1)read(333)bpara,bparasave
     if(cs_method==1)read(333)sapara_na,apara_na,phiind
     close(333)
#endif
     irest=irest+1 !!XY rolling rstart
     if(mype==0)write(gtcout,*)restart_fname,'read over'

  endif

end subroutine restart_native
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!----------------------------------------------------------------------------!
#ifdef ADIOS2
subroutine restart_adios2_write(restart_fname)
  use system_env
  use global_parameters
  use particle_array
  use field_array
  use conservative_scheme,only:sapara_na,apara_na,phiind
  !use simple_error_handler
  use mpi
  use adios2_common
  implicit none

  character(len=*), intent(in) :: restart_fname
  integer   :: ierr

  type(adios2_engine),   save :: adios2_restart_engine
  type(adios2_variable), save :: var_num_pes
  type(adios2_variable), save :: var_mi
  type(adios2_variable), save :: var_me
  type(adios2_variable), save :: var_rdtemi
  type(adios2_variable), save :: var_rddeni
  type(adios2_variable), save :: var_pfluxi
  type(adios2_variable), save :: var_phi
  type(adios2_variable), save :: var_phip00
  type(adios2_variable), save :: var_phi00
  type(adios2_variable), save :: var_zonali
  type(adios2_variable), save :: var_zonalci
  type(adios2_variable), save :: var_meshti
  type(adios2_variable), save :: var_meshte
  type(adios2_variable), save :: var_meshni
  type(adios2_variable), save :: var_meshne
  type(adios2_variable), save :: var_pmarki
  type(adios2_variable), save :: var_markeri
  type(adios2_variable), save :: var_markerit
  type(adios2_variable), save :: var_zion
  type(adios2_variable), save :: var_tfracn
  type(adios2_variable), save :: var_rdteme
  type(adios2_variable), save :: var_rddene
  type(adios2_variable), save :: var_pfluxe
  type(adios2_variable), save :: var_zonale
  type(adios2_variable), save :: var_zonalce
  type(adios2_variable), save :: var_pmarke
  type(adios2_variable), save :: var_markere
  type(adios2_variable), save :: var_tfrac
  type(adios2_variable), save :: var_zelectron
  type(adios2_variable), save :: var_densitye
  type(adios2_variable), save :: var_pressureepara
  type(adios2_variable), save :: var_pressureeperp
  type(adios2_variable), save :: var_phisave
  type(adios2_variable), save :: var_dnesave
  type(adios2_variable), save :: var_mf
  type(adios2_variable), save :: var_rdtemf
  type(adios2_variable), save :: var_rddenf
  type(adios2_variable), save :: var_zonalf
  type(adios2_variable), save :: var_zonalcf
  type(adios2_variable), save :: var_meshtf
  type(adios2_variable), save :: var_meshnf
  type(adios2_variable), save :: var_pfluxf
  type(adios2_variable), save :: var_pmarkf
  type(adios2_variable), save :: var_markerf
  type(adios2_variable), save :: var_zfast
  type(adios2_variable), save :: var_mfe
  type(adios2_variable), save :: var_fetfracn
  type(adios2_variable), save :: var_rdtemfe
  type(adios2_variable), save :: var_rddenfe
  type(adios2_variable), save :: var_zonalfe
  type(adios2_variable), save :: var_zonalcfe
  type(adios2_variable), save :: var_meshtfe
  type(adios2_variable), save :: var_meshnfe
  type(adios2_variable), save :: var_pfluxfe
  type(adios2_variable), save :: var_pmarkfe
  type(adios2_variable), save :: var_markerfe
  type(adios2_variable), save :: var_fetfrac
  type(adios2_variable), save :: var_zfaste
  type(adios2_variable), save :: var_apara00
  type(adios2_variable), save :: var_fluidne00
  type(adios2_variable), save :: var_apara
  type(adios2_variable), save :: var_sapara
  type(adios2_variable), save :: var_sdelapara
  type(adios2_variable), save :: var_fluidne
  type(adios2_variable), save :: var_sfluidne
  type(adios2_variable), save :: var_fluidue
  type(adios2_variable), save :: var_deltapsi
  type(adios2_variable), save :: var_sdeltapsi
  type(adios2_variable), save :: var_aparanl
  type(adios2_variable), save :: var_gradgaugef
  type(adios2_variable), save :: var_bpara
  type(adios2_variable), save :: var_bparasave
  type(adios2_variable), save :: var_sapara_na
  type(adios2_variable), save :: var_apara_na
  type(adios2_variable), save :: var_phiind

  integer   :: shape2(2), shape3(3)

  !! Initialize adios2 io object and declare variables
  if (.not.adios2_restart_io_initialized) then
    call adios2_declare_io (adios2_restart_io, adios2_obj, adios2_restart_io_id, ierr)
    call error_check(ierr, "adios2 restart_write declare_io", fatal=.true.)

    call adios2_define_variable (var_num_pes, adios2_restart_io, "num_pes", adios2_type_integer4, ierr)

    ! Create a 1D global array of mi 
    call adios2_define_variable (var_mi, adios2_restart_io, "mi", adios2_type_integer4 , 1,  &
                                 (/ 1_8 * numberpe /),                &  ! global shape
                                 (/ 1_8 * mype /),                &  ! offset
                                 (/ 1_8 /),                           &  ! count
                                 .true., ierr)
    call adios2_define_variable (var_me, adios2_restart_io, "me", adios2_type_integer4 , 1,  &
                                 (/ 1_8 * numberpe /),                &  ! global shape
                                 (/ 1_8 * mype /),                &  ! offset
                                 (/ 1_8 /),                           &  ! count
                                 .true., ierr)
    call adios2_define_variable (var_rdtemi, adios2_restart_io, "rdtemi", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(rdtemi), numberpe /),              &
                                 1_8 * (/ 0, mype /),              &
                                 1_8 * (/ shape(rdtemi), 1 /),                        &
                                 .true., ierr)
    call adios2_define_variable (var_rddeni, adios2_restart_io, "rddeni", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(rddeni), numberpe /),&
                                 1_8 * (/ 0, mype /),              &
                                 1_8 * (/ shape(rddeni), 1 /),                        &
                                 .true., ierr)
    call adios2_define_variable (var_pfluxi, adios2_restart_io, "pfluxi", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(pfluxi), numberpe /),              &
                                 1_8 * (/ 0, mype /),              &
                                 1_8 * (/ shape(pfluxi), 1 /),                        &
                                 .true., ierr)
    ! 2D arrays are combined by combining their large dimension
    shape2 = shape(phi)
    call adios2_define_variable (var_phi, adios2_restart_io, "phi", adios2_type_real_lk , 3,        &
                                 1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                 1_8 * (/ 0, 0, mype /), &
                                 1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
    call adios2_define_variable (var_phip00, adios2_restart_io, "phip00", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(phip00), numberpe /),              &
                                 1_8 * (/ 0, mype /),              &
                                 1_8 * (/ shape(phip00), 1 /),                        &
                                 .true., ierr)
    call adios2_define_variable (var_phi00, adios2_restart_io, "phi00", adios2_type_real_lk , 2,   &
                                 1_8 * (/ shape(phi00), numberpe /),              &
                                 1_8 * (/ 0, mype /),              &
                                 1_8 * (/ shape(phi00), 1 /),                         &
                                 .true., ierr)
    call adios2_define_variable (var_zonali, adios2_restart_io, "zonali", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(zonali), numberpe /),              &
                                 1_8 * (/ 0, mype /),              &
                                 1_8 * (/ shape(zonali), 1 /),                         &
                                 .true., ierr)
    call adios2_define_variable (var_zonalci, adios2_restart_io, "zonalci", adios2_type_real_lk , 2, &
                                 1_8 * (/ shape(zonalci), numberpe /),              &
                                 1_8 * (/ 0, mype /),             &
                                 1_8 * (/ shape(zonalci), 1 /),                         &
                                 .true., ierr)
    call adios2_define_variable (var_meshti, adios2_restart_io, "meshti", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(meshti), numberpe /),              &
                                 1_8 * (/ 0, mype /),             &
                                 1_8 * (/ shape(meshti), 1 /),                        &
                                 .true., ierr)
    call adios2_define_variable (var_meshte, adios2_restart_io, "meshte", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(meshte), numberpe /),              &
                                 1_8 * (/ 0, mype /),             &
                                 1_8 * (/ shape(meshte), 1 /),                         &
                                 .true., ierr)
    call adios2_define_variable (var_meshni, adios2_restart_io, "meshni", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(meshni), numberpe /),              &
                                 1_8 * (/ 0, mype /),             &
                                 1_8 * (/ shape(meshni), 1 /),                         &
                                 .true., ierr)
    call adios2_define_variable (var_meshne, adios2_restart_io, "meshne", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(meshne), numberpe /),              &
                                 1_8 * (/ 0, mype /),             &
                                 1_8 * (/ shape(meshne), 1 /),                         &
                                 .true., ierr)
    call adios2_define_variable (var_pmarki, adios2_restart_io, "pmarki", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(pmarki), numberpe /),              &
                                 1_8 * (/ 0, mype /),             &
                                 1_8 * (/ shape(pmarki), 1 /),                         &
                                 .true., ierr)
    call adios2_define_variable (var_markeri, adios2_restart_io, "markeri", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(markeri), numberpe /),              &
                                 1_8 * (/ 0, mype /),             &
                                 1_8 * (/ shape(markeri), 1 /),                         &
                                 .true., ierr)
    call adios2_define_variable (var_markerit, adios2_restart_io, "markerit", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(markerit), numberpe /),              &
                                 1_8 * (/ 0, mype /),             &
                                 1_8 * (/ shape(markerit), 1 /),                         &
                                 .true., ierr)
    ! Lets write the full zion instead of zion(:,1:mi)
    shape2 = shape(zion)
    call adios2_define_variable (var_zion, adios2_restart_io, "zion", adios2_type_real_lk , 3,        &
                                 1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                 1_8 * (/ 0, 0, mype /), &
                                 1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)

    if (nhybrid>0) then
      call adios2_define_variable (var_tfracn, adios2_restart_io, "tfracn", adios2_type_real_lk , 1,  &
                                   (/ 1_8 * numberpe /),                &  ! global shape
                                   (/ 1_8 * mype /),                &  ! offset
                                   (/ 1_8 /),                           &  ! count
                                   .true., ierr)
      call adios2_define_variable (var_rdteme, adios2_restart_io, "rdteme", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(rdteme), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(rdteme), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_rddene, adios2_restart_io, "rddene", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(rddene), numberpe /),&
                                 1_8 * (/ 0, mype /),              &
                                 1_8 * (/ shape(rddene), 1 /),&
                                 .true., ierr)
      call adios2_define_variable (var_pfluxe, adios2_restart_io, "pfluxe", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(pfluxe), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(pfluxe), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_zonale, adios2_restart_io, "zonale", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(zonale), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(zonale), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_zonalce, adios2_restart_io, "zonalce", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(zonalce), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(zonalce), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_pmarke, adios2_restart_io, "pmarke", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(pmarke), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(pmarke), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_markere, adios2_restart_io, "markere", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(markere), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(markere), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_tfrac, adios2_restart_io, "tfrac", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(tfrac), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(tfrac), 1 /),                         &
                                   .true., ierr)
      shape2 = shape(zelectron)
      call adios2_define_variable (var_zelectron, adios2_restart_io, "zelectron", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape2 = shape(densitye)
      call adios2_define_variable (var_densitye, adios2_restart_io, "densitye", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape2 = shape(pressureepara)
      call adios2_define_variable (var_pressureepara, adios2_restart_io, "pressureepara", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape2 = shape(pressureeperp)
      call adios2_define_variable (var_pressureeperp, adios2_restart_io, "pressureeperp", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape3 = shape(phisave)
      call adios2_define_variable (var_phisave, adios2_restart_io, "phisave", adios2_type_real_lk , 4,        &
                                   1_8 * (/ shape3(1), shape3(2), shape3(3), numberpe /), &
                                   1_8 * (/ 0, 0, 0, mype /), &
                                   1_8 * (/ shape3(1), shape3(2), shape3(3), 1 /), .true., ierr)
      shape3 = shape(dnesave)
      call adios2_define_variable (var_dnesave, adios2_restart_io, "dnesave", adios2_type_real_lk , 4,        &
                                   1_8 * (/ shape3(1), shape3(2), shape3(3), numberpe /), &
                                   1_8 * (/ 0, 0, 0, mype /), &
                                   1_8 * (/ shape3(1), shape3(2), shape3(3), 1 /), .true., ierr)
    endif
    
    if(fload>0) then
      call adios2_define_variable (var_mf, adios2_restart_io, "mf", adios2_type_integer4 , 1,  &
                                   (/ 1_8 * numberpe /),                &  ! global shape
                                   (/ 1_8 * mype /),                &  ! offset
                                   (/ 1_8 /),                           &  ! count
                                   .true., ierr)
      call adios2_define_variable (var_rdtemf, adios2_restart_io, "rdtemf", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(rdtemf), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(rdtemf), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_rddenf, adios2_restart_io, "rddenf", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(rddenf), numberpe /),&
                                 1_8 * (/ 0, mype /),              &
                                 1_8 * (/ shape(rddenf), 1 /),&
                                 .true., ierr)
      call adios2_define_variable (var_zonalf, adios2_restart_io, "zonalf", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(zonalf), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(zonalf), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_zonalcf, adios2_restart_io, "zonalcf", adios2_type_real_lk , 2, &
                                   1_8 * (/ shape(zonalcf), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(zonalcf), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_meshtf, adios2_restart_io, "meshtf", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(meshtf), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(meshtf), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_meshnf, adios2_restart_io, "meshnf", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(meshnf), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(meshnf), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_pfluxf, adios2_restart_io, "pfluxf", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(pfluxf), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(pfluxf), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_pmarkf, adios2_restart_io, "pmarkf", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(pmarkf), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(pmarkf), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_markerf, adios2_restart_io, "markerf", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(markerf), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(markerf), 1 /),                         &
                                   .true., ierr)
      ! Lets write the full zfast
      shape2 = shape(zfast)
      call adios2_define_variable (var_zfast, adios2_restart_io, "zfast", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
    endif
    
    if(feload>0) then
      call adios2_define_variable (var_mfe, adios2_restart_io, "mfe", adios2_type_integer4 , 1,  &
                                   (/ 1_8 * numberpe /),                &  ! global shape
                                   (/ 1_8 * mype /),                &  ! offset
                                   (/ 1_8 /),                           &  ! count
                                   .true., ierr)
      call adios2_define_variable (var_fetfracn, adios2_restart_io, "fetfracn", adios2_type_real_lk , 1,  &
                                   (/ 1_8 * numberpe /),                &  ! global shape
                                   (/ 1_8 * mype /),                &  ! offset
                                   (/ 1_8 /),                           &  ! count
                                   .true., ierr)
      call adios2_define_variable (var_rdtemfe, adios2_restart_io, "rdtemfe", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(rdtemfe), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(rdtemfe), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_rddenfe, adios2_restart_io, "rddenfe", adios2_type_real_lk , 2,  &
                                 1_8 * (/ shape(rddenfe), numberpe /),&
                                 1_8 * (/ 0, mype /),              &
                                 1_8 * (/ shape(rddenfe), 1 /),&
                                 .true., ierr)
      call adios2_define_variable (var_zonalfe, adios2_restart_io, "zonalfe", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(zonalfe), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(zonalfe), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_zonalcfe, adios2_restart_io, "zonalcfe", adios2_type_real_lk , 2, &
                                   1_8 * (/ shape(zonalcfe), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(zonalcfe), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_meshtfe, adios2_restart_io, "meshtfe", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(meshtfe), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(meshtfe), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_meshnfe, adios2_restart_io, "meshnfe", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(meshnfe), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(meshnfe), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_pfluxfe, adios2_restart_io, "pfluxfe", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(pfluxfe), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(pfluxfe), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_pmarkfe, adios2_restart_io, "pmarkfe", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(pmarkfe), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(pmarkfe), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_markerfe, adios2_restart_io, "markerfe", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(markerfe), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(markerfe), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_fetfrac, adios2_restart_io, "fetfrac", adios2_type_real_lk , 2,  &
                                   1_8 * (/ shape(fetfrac), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(fetfrac), 1 /),                         &
                                   .true., ierr)
      ! Lets write the full zfaste
      shape2 = shape(zfaste)
      call adios2_define_variable (var_zfaste, adios2_restart_io, "zfaste", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
    endif
    
    if(magnetic>0) then
      ! 1D - apara00, fluidne00,
      call adios2_define_variable (var_apara00, adios2_restart_io, "apara00", adios2_type_real_lk, 2,  &
                                   1_8 * (/ shape(apara00), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(apara00), 1 /),                         &
                                   .true., ierr)
      call adios2_define_variable (var_fluidne00, adios2_restart_io, "fluidne00", adios2_type_real_lk, 2,  &
                                   1_8 * (/ shape(fluidne00), numberpe /),              &
                                   1_8 * (/ 0, mype /),             &
                                   1_8 * (/ shape(fluidne00), 1 /),                         &
                                   .true., ierr)
      ! 2D apara,sapara,sdelapara,fluidne,sfluidne,fluidue,deltapsi,sdeltapsi,aparanl,
      shape2 = shape(apara)
      call adios2_define_variable (var_apara, adios2_restart_io, "apara", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape2 = shape(sapara)
      call adios2_define_variable (var_sapara, adios2_restart_io, "sapara", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape2 = shape(sdelapara)
      call adios2_define_variable (var_sdelapara, adios2_restart_io, "sdelapara", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape2 = shape(fluidne)
      call adios2_define_variable (var_fluidne, adios2_restart_io, "fluidne", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape2 = shape(sfluidne)
      call adios2_define_variable (var_sfluidne, adios2_restart_io, "sfluidne", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape2 = shape(fluidue)
      call adios2_define_variable (var_fluidue, adios2_restart_io, "fluidue", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape2 = shape(deltapsi)
      call adios2_define_variable (var_deltapsi, adios2_restart_io, "deltapsi", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape2 = shape(sdeltapsi)
      call adios2_define_variable (var_sdeltapsi, adios2_restart_io, "sdeltapsi", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape2 = shape(aparanl)
      call adios2_define_variable (var_aparanl, adios2_restart_io, "aparanl", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      ! 3D - gradgaugef
      shape3 = shape(gradgaugef)
      call adios2_define_variable (var_gradgaugef, adios2_restart_io, "gradgaugef", adios2_type_real_lk , 4,        &
                                   1_8 * (/ shape3(1), shape3(2), shape3(3), numberpe /), &
                                   1_8 * (/ 0, 0, 0, mype /), &
                                   1_8 * (/ shape3(1), shape3(2), shape3(3), 1 /), .true., ierr)
    endif

    if(deltabpara==1) then
      !bpara,bparasave
      shape2 = shape(bpara)
      call adios2_define_variable (var_bpara, adios2_restart_io, "bpara", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape3 = shape(bparasave)
      call adios2_define_variable (var_bparasave, adios2_restart_io, "bparasave", adios2_type_real_lk , 4,        &
                                   1_8 * (/ shape3(1), shape3(2), shape3(3), numberpe /), &
                                   1_8 * (/ 0, 0, 0, mype /), &
                                   1_8 * (/ shape3(1), shape3(2), shape3(3), 1 /), .true., ierr)
    endif

    if(cs_method==1) then
      !sapara_na,apara_na,phiind
      shape2 = shape(sapara_na)
      call adios2_define_variable (var_sapara_na, adios2_restart_io, "sapara_na", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape2 = shape(apara_na)
      call adios2_define_variable (var_apara_na, adios2_restart_io, "apara_na", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
      shape2 = shape(phiind)
      call adios2_define_variable (var_phiind, adios2_restart_io, "phiind", adios2_type_real_lk , 3,        &
                                   1_8 * (/ shape2(1), shape2(2), numberpe /), &
                                   1_8 * (/ 0, 0, mype /), &
                                   1_8 * (/ shape2(1), shape2(2), 1 /), .true., ierr)
    endif
    
    adios2_restart_io_initialized = .true.
  endif

  !! Open adios file. This returns an engine handle
  call adios2_open (adios2_restart_engine, adios2_restart_io, restart_fname, adios2_mode_write, mpi_comm_world, ierr)
  call error_check (ierr, "adios2_open in restart_write", fatal=.true.)

  call adios2_begin_step (adios2_restart_engine, adios2_step_mode_append, ierr)

  if (mype.eq.0) call adios2_put (adios2_restart_engine, var_num_pes, numberpe, ierr)

  !$acc update host(zion,zonali,zonalci)

  !! Write all variables
  call adios2_put (adios2_restart_engine, var_mi, mi, ierr)
  call adios2_put (adios2_restart_engine, var_me, me, ierr)
  call adios2_put (adios2_restart_engine, var_rdtemi, rdtemi, ierr)
  call adios2_put (adios2_restart_engine, var_rddeni, rddeni, ierr)
  call adios2_put (adios2_restart_engine, var_pfluxi, pfluxi, ierr)
  call adios2_put (adios2_restart_engine, var_phi, phi, ierr)
  call adios2_put (adios2_restart_engine, var_phip00, phip00, ierr)
  call adios2_put (adios2_restart_engine, var_phi00, phi00, ierr)
  call adios2_put (adios2_restart_engine, var_zonali, zonali, ierr)
  call adios2_put (adios2_restart_engine, var_zonalci, zonalci, ierr)
  call adios2_put (adios2_restart_engine, var_meshti, meshti, ierr)
  call adios2_put (adios2_restart_engine, var_meshte, meshte, ierr)
  call adios2_put (adios2_restart_engine, var_meshni, meshni, ierr)
  call adios2_put (adios2_restart_engine, var_meshne, meshne, ierr)
  call adios2_put (adios2_restart_engine, var_pmarki, pmarki, ierr)
  call adios2_put (adios2_restart_engine, var_markeri, markeri, ierr)
  call adios2_put (adios2_restart_engine, var_markerit, markerit, ierr)
  call adios2_put (adios2_restart_engine, var_zion, zion, ierr)

  if (nhybrid>0) then
    !$acc update host(zelectron,zonale,zonalce)
    call adios2_put (adios2_restart_engine, var_tfracn, tfracn, ierr)
    call adios2_put (adios2_restart_engine, var_rdteme, rdteme, ierr)
    call adios2_put (adios2_restart_engine, var_rddene, rddene, ierr)
    call adios2_put (adios2_restart_engine, var_pfluxe, pfluxe, ierr)
    call adios2_put (adios2_restart_engine, var_zonale, zonale, ierr)
    call adios2_put (adios2_restart_engine, var_zonalce, zonalce, ierr)
    call adios2_put (adios2_restart_engine, var_pmarke, pmarke, ierr)
    call adios2_put (adios2_restart_engine, var_markere, markere, ierr)
    call adios2_put (adios2_restart_engine, var_tfrac, tfrac, ierr)
    call adios2_put (adios2_restart_engine, var_zelectron, zelectron, ierr)
    call adios2_put (adios2_restart_engine, var_densitye, densitye, ierr)
    call adios2_put (adios2_restart_engine, var_pressureepara, pressureepara, ierr)
    call adios2_put (adios2_restart_engine, var_pressureeperp, pressureeperp, ierr)
    call adios2_put (adios2_restart_engine, var_phisave, phisave, ierr)
    call adios2_put (adios2_restart_engine, var_dnesave, dnesave, ierr)
  endif
  
  if(fload>0) then
    !$acc update host(zfast,zonalf,zonalcf)
    call adios2_put (adios2_restart_engine, var_mf, mf, ierr)
    call adios2_put (adios2_restart_engine, var_rdtemf, rdtemf, ierr)
    call adios2_put (adios2_restart_engine, var_rddenf, rddenf, ierr)
    call adios2_put (adios2_restart_engine, var_zonalf, zonalf, ierr)
    call adios2_put (adios2_restart_engine, var_zonalcf, zonalcf, ierr)
    call adios2_put (adios2_restart_engine, var_meshtf, meshtf, ierr)
    call adios2_put (adios2_restart_engine, var_meshnf, meshnf, ierr)
    call adios2_put (adios2_restart_engine, var_pfluxf, pfluxf, ierr)
    call adios2_put (adios2_restart_engine, var_pmarkf, pmarkf, ierr)
    call adios2_put (adios2_restart_engine, var_markerf, markerf, ierr)
    call adios2_put (adios2_restart_engine, var_zfast, zfast, ierr)
  endif
  
  if(feload>0) then
    !$acc update host(zfaste,zonalfe,zonalcfe)
    call adios2_put (adios2_restart_engine, var_mfe, mfe, ierr)
    call adios2_put (adios2_restart_engine, var_fetfracn, fetfracn, ierr)
    call adios2_put (adios2_restart_engine, var_rdtemfe, rdtemfe, ierr)
    call adios2_put (adios2_restart_engine, var_rddenfe, rddenfe, ierr)
    call adios2_put (adios2_restart_engine, var_zonalfe, zonalfe, ierr)
    call adios2_put (adios2_restart_engine, var_zonalcfe, zonalcfe, ierr)
    call adios2_put (adios2_restart_engine, var_meshtfe, meshtfe, ierr)
    call adios2_put (adios2_restart_engine, var_meshnfe, meshnfe, ierr)
    call adios2_put (adios2_restart_engine, var_pfluxfe, pfluxfe, ierr)
    call adios2_put (adios2_restart_engine, var_pmarkfe, pmarkfe, ierr)
    call adios2_put (adios2_restart_engine, var_markerfe, markerfe, ierr)
    call adios2_put (adios2_restart_engine, var_fetfrac, fetfrac, ierr)
    call adios2_put (adios2_restart_engine, var_zfaste, zfaste, ierr)
  endif
  
  if(magnetic>0) then
    call adios2_put (adios2_restart_engine, var_apara00, apara00, ierr)
    call adios2_put (adios2_restart_engine, var_fluidne00, fluidne00, ierr)
    call adios2_put (adios2_restart_engine, var_apara, apara, ierr)
    call adios2_put (adios2_restart_engine, var_sapara, sapara, ierr)
    call adios2_put (adios2_restart_engine, var_sdelapara, sdelapara, ierr)
    call adios2_put (adios2_restart_engine, var_fluidne, fluidne, ierr)
    call adios2_put (adios2_restart_engine, var_sfluidne, sfluidne, ierr)
    call adios2_put (adios2_restart_engine, var_fluidue, fluidue, ierr)
    call adios2_put (adios2_restart_engine, var_deltapsi, deltapsi, ierr)
    call adios2_put (adios2_restart_engine, var_sdeltapsi, sdeltapsi, ierr)
    call adios2_put (adios2_restart_engine, var_aparanl, aparanl, ierr)
    call adios2_put (adios2_restart_engine, var_gradgaugef, gradgaugef, ierr)
  endif

  if(deltabpara==1) then
    call adios2_put (adios2_restart_engine, var_bpara, bpara, ierr)
    call adios2_put (adios2_restart_engine, var_bparasave, bparasave, ierr)
  endif

  if(cs_method==1) then
    call adios2_put (adios2_restart_engine, var_sapara_na, sapara_na, ierr)
    call adios2_put (adios2_restart_engine, var_apara_na, apara_na, ierr)
    call adios2_put (adios2_restart_engine, var_phiind, phiind, ierr)
  endif
    
  call adios2_end_step (adios2_restart_engine, ierr)

  !! Close file object
  call adios2_close (adios2_restart_engine, ierr)
  call error_check (ierr, "adios2_close in restart_write", fatal=.true.)
end subroutine restart_adios2_write


!----------------------------------------------------------------------------!
subroutine restart_adios2_read(restart_fname)
  use system_env
  use global_parameters
  use particle_array
  use field_array
  use conservative_scheme,only:sapara_na,apara_na,phiind
  use simple_error_handler
  use mpi
  use adios2_common
  implicit none

  character(len=*), intent(in)  :: restart_fname
  integer   :: step_status, ierr
  integer   :: shape2(2), shape3(3)
  
  type(adios2_io)       :: adios2_restart_read_io
  type(adios2_engine)   :: adios2_restart_engine
  type(adios2_variable) :: var_mi
  type(adios2_variable) :: var_me
  type(adios2_variable) :: var_rdtemi
  type(adios2_variable) :: var_rddeni
  type(adios2_variable) :: var_pfluxi
  type(adios2_variable) :: var_phi
  type(adios2_variable) :: var_phip00
  type(adios2_variable) :: var_phi00
  type(adios2_variable) :: var_zonali
  type(adios2_variable) :: var_zonalci
  type(adios2_variable) :: var_meshti
  type(adios2_variable) :: var_meshte
  type(adios2_variable) :: var_meshni
  type(adios2_variable) :: var_meshne
  type(adios2_variable) :: var_pmarki
  type(adios2_variable) :: var_markeri
  type(adios2_variable) :: var_markerit
  type(adios2_variable) :: var_zion
  type(adios2_variable) :: var_tfracn
  type(adios2_variable) :: var_rdteme
  type(adios2_variable) :: var_rddene
  type(adios2_variable) :: var_pfluxe
  type(adios2_variable) :: var_zonale
  type(adios2_variable) :: var_zonalce
  type(adios2_variable) :: var_pmarke
  type(adios2_variable) :: var_markere
  type(adios2_variable) :: var_tfrac
  type(adios2_variable) :: var_zelectron
  type(adios2_variable) :: var_densitye
  type(adios2_variable) :: var_pressureepara
  type(adios2_variable) :: var_pressureeperp
  type(adios2_variable) :: var_phisave
  type(adios2_variable) :: var_dnesave
  type(adios2_variable) :: var_mf
  type(adios2_variable) :: var_rdtemf
  type(adios2_variable) :: var_rddenf
  type(adios2_variable) :: var_zonalf
  type(adios2_variable) :: var_zonalcf
  type(adios2_variable) :: var_meshtf
  type(adios2_variable) :: var_meshnf
  type(adios2_variable) :: var_pfluxf
  type(adios2_variable) :: var_pmarkf
  type(adios2_variable) :: var_markerf
  type(adios2_variable) :: var_zfast
  type(adios2_variable) :: var_mfe
  type(adios2_variable) :: var_fetfracn
  type(adios2_variable) :: var_rdtemfe
  type(adios2_variable) :: var_rddenfe
  type(adios2_variable) :: var_zonalfe
  type(adios2_variable) :: var_zonalcfe
  type(adios2_variable) :: var_meshtfe
  type(adios2_variable) :: var_meshnfe
  type(adios2_variable) :: var_pfluxfe
  type(adios2_variable) :: var_pmarkfe
  type(adios2_variable) :: var_markerfe
  type(adios2_variable) :: var_fetfrac
  type(adios2_variable) :: var_zfaste
  type(adios2_variable) :: var_apara00
  type(adios2_variable) :: var_fluidne00
  type(adios2_variable) :: var_apara
  type(adios2_variable) :: var_sapara
  type(adios2_variable) :: var_sdelapara
  type(adios2_variable) :: var_fluidne
  type(adios2_variable) :: var_sfluidne
  type(adios2_variable) :: var_fluidue
  type(adios2_variable) :: var_deltapsi
  type(adios2_variable) :: var_sdeltapsi
  type(adios2_variable) :: var_aparanl
  type(adios2_variable) :: var_gradgaugef
  type(adios2_variable) :: var_bpara
  type(adios2_variable) :: var_bparasave
  type(adios2_variable) :: var_sapara_na
  type(adios2_variable) :: var_apara_na
  type(adios2_variable) :: var_phiind

  call mpi_comm_size (mpi_comm_world, numberpe, ierr)
  call mpi_comm_rank (mpi_comm_world, mype, ierr)

  !! Initialize adios2 io object and declare variables
  call adios2_declare_io (adios2_restart_read_io, adios2_obj, adios2_restart_io_read_id, ierr)
  call error_check(ierr, "adios2 restart_read declare_io", fatal=.true.)

  !! Open adios file. This returns an engine handle
  call adios2_open (adios2_restart_engine, adios2_restart_read_io, restart_fname, adios2_mode_read, mpi_comm_world, ierr)
  call error_check (ierr, "adios2_open in restart_read", fatal=.true.)

  call adios2_begin_step(adios2_restart_engine, adios2_step_mode_read, 0.0, step_status, ierr)

  call adios2_inquire_variable (var_mi, adios2_restart_read_io, "mi", ierr)
  call adios2_set_selection (var_mi, 1, (/ 1_8 * mype /), (/ 1_8 /), ierr)

  call adios2_inquire_variable (var_me, adios2_restart_read_io, "me", ierr)
  call adios2_set_selection (var_me, 1, (/ 1_8 * mype /), (/ 1_8 /), ierr)

  call adios2_inquire_variable (var_rdtemi, adios2_restart_read_io, "rdtemi", ierr)
  call adios2_set_selection (var_rdtemi, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(rdtemi), 1 /), ierr)

  call adios2_inquire_variable (var_rddeni, adios2_restart_read_io, "rddeni", ierr)
  call adios2_set_selection (var_rddeni, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(rddeni), 1 /), ierr)
  
  call adios2_inquire_variable (var_pfluxi, adios2_restart_read_io, "pfluxi", ierr)
  call adios2_set_selection (var_pfluxi, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(pfluxi), 1 /), ierr)
  
  shape2 = shape(phi)
  call adios2_inquire_variable (var_phi, adios2_restart_read_io, "phi", ierr)
  call adios2_set_selection (var_phi, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1 /), ierr)
  
  call adios2_inquire_variable (var_phip00, adios2_restart_read_io, "phip00", ierr)
  call adios2_set_selection (var_phip00, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(phip00), 1 /), ierr)
  
  call adios2_inquire_variable (var_phi00, adios2_restart_read_io, "phi00", ierr)
  call adios2_set_selection (var_phi00, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(phi00), 1/), ierr)
  
  call adios2_inquire_variable (var_zonali, adios2_restart_read_io, "zonali", ierr)
  call adios2_set_selection (var_zonali, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(zonali), 1/), ierr)
  
  call adios2_inquire_variable (var_zonalci, adios2_restart_read_io, "zonalci", ierr)
  call adios2_set_selection (var_zonalci, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(zonalci), 1/), ierr)
  
  call adios2_inquire_variable (var_meshti, adios2_restart_read_io, "meshti", ierr)
  call adios2_set_selection (var_meshti, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(meshti), 1/), ierr)
  
  call adios2_inquire_variable (var_meshte, adios2_restart_read_io, "meshte", ierr)
  call adios2_set_selection (var_meshte, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(meshte), 1/), ierr)
  
  call adios2_inquire_variable (var_meshni, adios2_restart_read_io, "meshni", ierr)
  call adios2_set_selection (var_meshni, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(meshni), 1/), ierr)
  
  call adios2_inquire_variable (var_meshne, adios2_restart_read_io, "meshne", ierr)
  call adios2_set_selection (var_meshne, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(meshne), 1/), ierr)
  
  call adios2_inquire_variable (var_pmarki, adios2_restart_read_io, "pmarki", ierr)
  call adios2_set_selection (var_pmarki, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(pmarki), 1/), ierr)
  
  call adios2_inquire_variable (var_markeri, adios2_restart_read_io, "markeri", ierr)
  call adios2_set_selection (var_markeri, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(markeri), 1/), ierr)
  
  call adios2_inquire_variable (var_markerit, adios2_restart_read_io, "markerit", ierr)
  call adios2_set_selection (var_markerit, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(markerit), 1/), ierr)
  
  shape2 = shape(zion)
  call adios2_inquire_variable (var_zion, adios2_restart_read_io, "zion", ierr)
  call adios2_set_selection (var_zion, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

  if (nhybrid>0) then
  
    call adios2_inquire_variable (var_tfracn, adios2_restart_read_io, "tfracn", ierr)
    call adios2_set_selection (var_tfracn, 1, (/ 1_8 * mype /), (/ 1_8 /), ierr)
    
    call adios2_inquire_variable (var_rdteme, adios2_restart_read_io, "rdteme", ierr)
    call adios2_set_selection (var_rdteme, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(rdteme), 1/), ierr)

    call adios2_inquire_variable (var_rddene, adios2_restart_read_io, "rddene", ierr)
    call adios2_set_selection (var_rddene, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(rddene), 1 /), ierr)
    
    call adios2_inquire_variable (var_pfluxe, adios2_restart_read_io, "pfluxe", ierr)
    call adios2_set_selection (var_pfluxe, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(pfluxe), 1/), ierr)
    
    call adios2_inquire_variable (var_zonale, adios2_restart_read_io, "zonale", ierr)
    call adios2_set_selection (var_zonale, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(zonale), 1/), ierr)
    
    call adios2_inquire_variable (var_zonalce, adios2_restart_read_io, "zonalce", ierr)
    call adios2_set_selection (var_zonalce, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(zonalce), 1/), ierr)
    
    call adios2_inquire_variable (var_pmarke, adios2_restart_read_io, "pmarke", ierr)
    call adios2_set_selection (var_pmarke, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(pmarke), 1/), ierr)
    
    call adios2_inquire_variable (var_markere, adios2_restart_read_io, "markere", ierr)
    call adios2_set_selection (var_markere, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(markere), 1/), ierr)
    
    call adios2_inquire_variable (var_tfrac, adios2_restart_read_io, "tfrac", ierr)
    call adios2_set_selection (var_tfrac, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(tfrac), 1/), ierr)
    
    shape2 = shape(zelectron)
    call adios2_inquire_variable (var_zelectron, adios2_restart_read_io, "zelectron", ierr)
    call adios2_set_selection (var_zelectron, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape2 = shape(densitye)
    call adios2_inquire_variable (var_densitye, adios2_restart_read_io, "densitye", ierr)
    call adios2_set_selection (var_densitye, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape2 = shape(pressureepara)
    call adios2_inquire_variable (var_pressureepara, adios2_restart_read_io, "pressureepara", ierr)
    call adios2_set_selection (var_pressureepara, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape2 = shape(pressureeperp)    
    call adios2_inquire_variable (var_pressureeperp, adios2_restart_read_io, "pressureeperp", ierr)
    call adios2_set_selection (var_pressureeperp, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape3 = shape(phisave)    
    call adios2_inquire_variable (var_phisave, adios2_restart_read_io, "phisave", ierr)
    call adios2_set_selection (var_phisave, 4, 1_8 * (/ 0, 0, 0, mype /), 1_8 * (/ shape3(1), shape3(2), shape3(3), 1/), ierr)

    shape3 = shape(dnesave)    
    call adios2_inquire_variable (var_dnesave, adios2_restart_read_io, "dnesave", ierr)
    call adios2_set_selection (var_dnesave, 4, 1_8 * (/ 0, 0, 0, mype /), 1_8 * (/ shape3(1), shape3(2), shape3(3), 1/), ierr)
  endif
  
  if(fload>0) then
    
    call adios2_inquire_variable (var_mf, adios2_restart_read_io, "mf", ierr)
    call adios2_set_selection (var_mf, 1, (/ 1_8 * mype /), (/ 1_8 /), ierr)
    
    call adios2_inquire_variable (var_rdtemf, adios2_restart_read_io, "rdtemf", ierr)
    call adios2_set_selection (var_rdtemf, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(rdtemf), 1/), ierr)

    call adios2_inquire_variable (var_rddenf, adios2_restart_read_io, "rddenf", ierr)
    call adios2_set_selection (var_rddenf, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(rddenf), 1 /), ierr)
    
    call adios2_inquire_variable (var_zonalf, adios2_restart_read_io, "zonalf", ierr)
    call adios2_set_selection (var_zonalf, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(zonalf), 1/), ierr)
    
    call adios2_inquire_variable (var_zonalcf, adios2_restart_read_io, "zonalcf", ierr)
    call adios2_set_selection (var_zonalcf, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(zonalcf), 1/), ierr)
    
    call adios2_inquire_variable (var_meshtf, adios2_restart_read_io, "meshtf", ierr)
    call adios2_set_selection (var_meshtf, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(meshtf), 1/), ierr)
    
    call adios2_inquire_variable (var_meshnf, adios2_restart_read_io, "meshnf", ierr)
    call adios2_set_selection (var_meshnf, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(meshnf), 1/), ierr)
    
    call adios2_inquire_variable (var_pfluxf, adios2_restart_read_io, "pfluxf", ierr)
    call adios2_set_selection (var_pfluxf, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(pfluxf), 1/), ierr)
    
    call adios2_inquire_variable (var_pmarkf, adios2_restart_read_io, "pmarkf", ierr)
    call adios2_set_selection (var_pmarkf, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(pmarkf), 1/), ierr)
    
    call adios2_inquire_variable (var_markerf, adios2_restart_read_io, "markerf", ierr)
    call adios2_set_selection (var_markerf, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(markerf), 1/), ierr)
    
    shape2 = shape(zfast)
    call adios2_inquire_variable (var_zfast, adios2_restart_read_io, "zfast", ierr)
    call adios2_set_selection (var_zfast, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)
  endif
  
  if(feload>0) then
    
    call adios2_inquire_variable (var_mfe, adios2_restart_read_io, "mfe", ierr)
    call adios2_set_selection (var_mfe, 1, (/ 1_8 * mype /), (/ 1_8 /), ierr)
    
    call adios2_inquire_variable (var_fetfracn, adios2_restart_read_io, "fetfracn", ierr)
    call adios2_set_selection (var_fetfracn, 1, (/ 1_8 * mype /), (/ 1_8 /), ierr)
    
    call adios2_inquire_variable (var_rdtemfe, adios2_restart_read_io, "rdtemfe", ierr)
    call adios2_set_selection (var_rdtemfe, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(rdtemfe), 1/), ierr)

    call adios2_inquire_variable (var_rddenfe, adios2_restart_read_io, "rddenfe", ierr)
    call adios2_set_selection (var_rddenfe, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(rddenfe), 1 /), ierr)
    
    call adios2_inquire_variable (var_zonalfe, adios2_restart_read_io, "zonalfe", ierr)
    call adios2_set_selection (var_zonalfe, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(zonalfe), 1/), ierr)
    
    call adios2_inquire_variable (var_zonalcfe, adios2_restart_read_io, "zonalcfe", ierr)
    call adios2_set_selection (var_zonalcfe, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(zonalcfe), 1/), ierr)
    
    call adios2_inquire_variable (var_meshtfe, adios2_restart_read_io, "meshtfe", ierr)
    call adios2_set_selection (var_meshtfe, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(meshtfe), 1/), ierr)
    
    call adios2_inquire_variable (var_meshnfe, adios2_restart_read_io, "meshnfe", ierr)
    call adios2_set_selection (var_meshnfe, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(meshnfe), 1/), ierr)
    
    call adios2_inquire_variable (var_pfluxfe, adios2_restart_read_io, "pfluxfe", ierr)
    call adios2_set_selection (var_pfluxfe, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(pfluxfe), 1/), ierr)
    
    call adios2_inquire_variable (var_pmarkfe, adios2_restart_read_io, "pmarkfe", ierr)
    call adios2_set_selection (var_pmarkfe, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(pmarkfe), 1/), ierr)
    
    call adios2_inquire_variable (var_markerfe, adios2_restart_read_io, "markerfe", ierr)
    call adios2_set_selection (var_markerfe, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(markerfe), 1/), ierr)
    
    call adios2_inquire_variable (var_fetfrac, adios2_restart_read_io, "fetfrac", ierr)
    call adios2_set_selection (var_fetfrac, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(fetfrac), 1/), ierr)
    
    shape2 = shape(zfaste)    
    call adios2_inquire_variable (var_zfaste, adios2_restart_read_io, "zfaste", ierr)
    call adios2_set_selection (var_zfaste, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)
  endif
  
  if(magnetic>0) then
    
    call adios2_inquire_variable (var_apara00, adios2_restart_read_io, "apara00", ierr)
    call adios2_set_selection (var_apara00, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(apara00), 1/), ierr)
    
    call adios2_inquire_variable (var_fluidne00, adios2_restart_read_io, "fluidne00", ierr)
    call adios2_set_selection (var_fluidne00, 2, 1_8 * (/ 0, mype /), 1_8 * (/ shape(fluidne00), 1/), ierr)
    
    shape2 = shape(apara)    
    call adios2_inquire_variable (var_apara, adios2_restart_read_io, "apara", ierr)
    call adios2_set_selection (var_apara, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape2 = shape(sapara)    
    call adios2_inquire_variable (var_sapara, adios2_restart_read_io, "sapara", ierr)
    call adios2_set_selection (var_sapara, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape2 = shape(sdelapara)    
    call adios2_inquire_variable (var_sdelapara, adios2_restart_read_io, "sdelapara", ierr)
    call adios2_set_selection (var_sdelapara, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape2 = shape(fluidne)    
    call adios2_inquire_variable (var_fluidne, adios2_restart_read_io, "fluidne", ierr)
    call adios2_set_selection (var_fluidne, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape2 = shape(sfluidne)    
    call adios2_inquire_variable (var_sfluidne, adios2_restart_read_io, "sfluidne", ierr)
    call adios2_set_selection (var_sfluidne, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape2 = shape(fluidue)    
    call adios2_inquire_variable (var_fluidue, adios2_restart_read_io, "fluidue", ierr)
    call adios2_set_selection (var_fluidue, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape2 = shape(deltapsi)    
    call adios2_inquire_variable (var_deltapsi, adios2_restart_read_io, "deltapsi", ierr)
    call adios2_set_selection (var_deltapsi, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape2 = shape(sdeltapsi)    
    call adios2_inquire_variable (var_sdeltapsi, adios2_restart_read_io, "sdeltapsi", ierr)
    call adios2_set_selection (var_sdeltapsi, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape2 = shape(aparanl)    
    call adios2_inquire_variable (var_aparanl, adios2_restart_read_io, "aparanl", ierr)
    call adios2_set_selection (var_aparanl, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape3 = shape(gradgaugef)    
    call adios2_inquire_variable (var_gradgaugef, adios2_restart_read_io, "gradgaugef", ierr)
    call adios2_set_selection (var_gradgaugef, 4, 1_8 * (/ 0, 0, 0, mype /), 1_8 * (/ shape3(1), shape3(2), shape3(3), 1/), ierr)
  endif
    
  if(deltabpara==1) then
    !bpara,bparasave
    shape2 = shape(bpara)
    call adios2_inquire_variable (var_bpara, adios2_restart_read_io, "bpara", ierr)
    call adios2_set_selection (var_bpara, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape3 = shape(bparasave)    
    call adios2_inquire_variable (var_bparasave, adios2_restart_read_io, "bparasave", ierr)
    call adios2_set_selection (var_bparasave, 4, 1_8 * (/ 0, 0, 0, mype /), 1_8 * (/ shape3(1), shape3(2), shape3(3), 1/), ierr)
  endif
    
  if(cs_method==1) then
    !sapara_na,apara_na,phiind
    shape2 = shape(sapara_na)    
    call adios2_inquire_variable (var_sapara_na, adios2_restart_read_io, "sapara_na", ierr)
    call adios2_set_selection (var_sapara_na, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape2 = shape(apara_na)    
    call adios2_inquire_variable (var_apara_na, adios2_restart_read_io, "apara_na", ierr)
    call adios2_set_selection (var_apara_na, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)

    shape2 = shape(phiind)    
    call adios2_inquire_variable (var_phiind, adios2_restart_read_io, "phiind", ierr)
    call adios2_set_selection (var_phiind, 3, 1_8 * (/ 0, 0, mype /), 1_8 * (/ shape2(1), shape2(2), 1/), ierr)
  endif

  !! Perform reads
  call adios2_get (adios2_restart_engine, var_mi, mi, ierr)
  call adios2_get (adios2_restart_engine, var_me, me, ierr)
  call adios2_get (adios2_restart_engine, var_rdtemi, rdtemi, ierr)
  call adios2_get (adios2_restart_engine, var_rddeni, rddeni, ierr)
  call adios2_get (adios2_restart_engine, var_pfluxi, pfluxi, ierr)
  call adios2_get (adios2_restart_engine, var_phi, phi, ierr)
  call adios2_get (adios2_restart_engine, var_phip00, phip00, ierr)
  call adios2_get (adios2_restart_engine, var_phi00, phi00, ierr)
  call adios2_get (adios2_restart_engine, var_zonali, zonali, ierr)
  call adios2_get (adios2_restart_engine, var_zonalci, zonalci, ierr)
  call adios2_get (adios2_restart_engine, var_meshti, meshti, ierr)
  call adios2_get (adios2_restart_engine, var_meshte, meshte, ierr)
  call adios2_get (adios2_restart_engine, var_meshni, meshni, ierr)
  call adios2_get (adios2_restart_engine, var_meshne, meshne, ierr)
  call adios2_get (adios2_restart_engine, var_pmarki, pmarki, ierr)
  call adios2_get (adios2_restart_engine, var_markeri, markeri, ierr)
  call adios2_get (adios2_restart_engine, var_markerit, markerit, ierr)
  call adios2_get (adios2_restart_engine, var_zion, zion, ierr)

  if (nhybrid>0) then
    call adios2_get (adios2_restart_engine, var_tfracn, tfracn, ierr)
    call adios2_get (adios2_restart_engine, var_rdteme, rdteme, ierr)
    call adios2_get (adios2_restart_engine, var_rddene, rddene, ierr)
    call adios2_get (adios2_restart_engine, var_pfluxe, pfluxe, ierr)
    call adios2_get (adios2_restart_engine, var_zonale, zonale, ierr)
    call adios2_get (adios2_restart_engine, var_zonalce, zonalce, ierr)
    call adios2_get (adios2_restart_engine, var_pmarke, pmarke, ierr)
    call adios2_get (adios2_restart_engine, var_markere, markere, ierr)
    call adios2_get (adios2_restart_engine, var_tfrac, tfrac, ierr)
    call adios2_get (adios2_restart_engine, var_zelectron, zelectron, ierr)
    call adios2_get (adios2_restart_engine, var_densitye, densitye, ierr)
    call adios2_get (adios2_restart_engine, var_pressureepara, pressureepara, ierr)
    call adios2_get (adios2_restart_engine, var_pressureeperp, pressureeperp, ierr)
    call adios2_get (adios2_restart_engine, var_phisave, phisave, ierr)
    call adios2_get (adios2_restart_engine, var_dnesave, dnesave, ierr)
    !$acc update device(phisave,dnesave,densitye,pressureepara,pressureeperp)
  endif

  if (fload > 0) then
    call adios2_get (adios2_restart_engine, var_mf, mf, ierr)
    call adios2_get (adios2_restart_engine, var_rdtemf, rdtemf, ierr)
    call adios2_get (adios2_restart_engine, var_rddenf, rddenf, ierr)
    call adios2_get (adios2_restart_engine, var_zonalf, zonalf, ierr)
    call adios2_get (adios2_restart_engine, var_zonalcf, zonalcf, ierr)
    call adios2_get (adios2_restart_engine, var_meshtf, meshtf, ierr)
    call adios2_get (adios2_restart_engine, var_meshnf, meshnf, ierr)
    call adios2_get (adios2_restart_engine, var_pfluxf, pfluxf, ierr)
    call adios2_get (adios2_restart_engine, var_pmarkf, pmarkf, ierr)
    call adios2_get (adios2_restart_engine, var_markerf, markerf, ierr)
    call adios2_get (adios2_restart_engine, var_zfast, zfast, ierr)
  endif

  if (feload > 0) then
    call adios2_get (adios2_restart_engine, var_mfe, mfe, ierr)
    call adios2_get (adios2_restart_engine, var_fetfracn, fetfracn, ierr)
    call adios2_get (adios2_restart_engine, var_rdtemfe, rdtemfe, ierr)
    call adios2_get (adios2_restart_engine, var_rddenfe, rddenfe, ierr)
    call adios2_get (adios2_restart_engine, var_zonalfe, zonalfe, ierr)
    call adios2_get (adios2_restart_engine, var_zonalcfe, zonalcfe, ierr)
    call adios2_get (adios2_restart_engine, var_meshtfe, meshtfe, ierr)
    call adios2_get (adios2_restart_engine, var_meshnfe, meshnfe, ierr)
    call adios2_get (adios2_restart_engine, var_pfluxfe, pfluxfe, ierr)
    call adios2_get (adios2_restart_engine, var_pmarkfe, pmarkfe, ierr)
    call adios2_get (adios2_restart_engine, var_markerfe, markerfe, ierr)
    call adios2_get (adios2_restart_engine, var_fetfrac, fetfrac, ierr)
    call adios2_get (adios2_restart_engine, var_zfaste, zfaste, ierr)
  endif

  if (magnetic > 0) then
    call adios2_get (adios2_restart_engine, var_apara00, apara00, ierr)
    call adios2_get (adios2_restart_engine, var_fluidne00, fluidne00, ierr)
    call adios2_get (adios2_restart_engine, var_apara, apara, ierr)
    call adios2_get (adios2_restart_engine, var_sapara, sapara, ierr)
    call adios2_get (adios2_restart_engine, var_sdelapara, sdelapara, ierr)
    call adios2_get (adios2_restart_engine, var_fluidne, fluidne, ierr)
    call adios2_get (adios2_restart_engine, var_sfluidne, sfluidne, ierr)
    call adios2_get (adios2_restart_engine, var_fluidue, fluidue, ierr)
    call adios2_get (adios2_restart_engine, var_deltapsi, deltapsi, ierr)
    call adios2_get (adios2_restart_engine, var_sdeltapsi, sdeltapsi, ierr)
    call adios2_get (adios2_restart_engine, var_aparanl, aparanl, ierr)
    call adios2_get (adios2_restart_engine, var_gradgaugef, gradgaugef, ierr)
  endif

  if(deltabpara==1) then
    call adios2_get (adios2_restart_engine, var_bpara, bpara, ierr)
    call adios2_get (adios2_restart_engine, var_bparasave, bparasave, ierr)
  endif

  if(cs_method==1) then
    call adios2_get (adios2_restart_engine, var_sapara_na, sapara_na, ierr)
    call adios2_get (adios2_restart_engine, var_apara_na, apara_na, ierr)
    call adios2_get (adios2_restart_engine, var_phiind, phiind, ierr)
  endif

  call adios2_end_step (adios2_restart_engine, ierr)

  !! Close file object
  call adios2_close (adios2_restart_engine, ierr)
end subroutine

!end conditional compilation for ADIOS2
#endif

!----------------------------------------------------------------------------!
subroutine restart_io(iop)
  use system_env
  use utility, only: check_allocation, mstat
  use global_parameters
  use particle_array
  use field_array
  use data_type
  use particle_tracking
  use conservative_scheme
  implicit none
  integer mrequest, mfmode
  integer i,j,k,m,subsize,startidx,endidx,ndata
  integer mquantity,mflx,n_mode,mstepfinal,noutputs
  real(lk) dum
  character(*),intent(in)::iop
  character(len=80) :: restart_fname
  character(len=18) cdum
  integer :: save_restart_files
  INTEGER(KIND=MPI_OFFSET_KIND) mype_filesize, sum_filesize

  if(iop/="read" .and. iop/="write")then
     write(gtcout,*)'*** subroutine restart_io (iop <> "read" or "write")',iop
     call GTC_ABORT
     return
  endif
  !$acc update host(apara,sapara,sdelapara,fluidne,fluidue,sfluidne,deltapsi,sdeltapsi)
  !$acc update host(aparanl,apara00,fluidne00,gradgaugef)
  !$acc update host(bpara,sapara_na,apara_na,phiind)


  call restart_native(iop)

     !$acc update device(apara,sapara,sdelapara,fluidne,sfluidne,fluidue,deltapsi,sdeltapsi)
     !$acc update device(aparanl,apara00,fluidne00,gradgaugef)
     !$acc update device(bpara,sapara_na,apara_na,phiind)
end subroutine restart_io
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine restart_hist
  use global_parameters
  use utility
  implicit none
  integer i,ii(6),mstepfinal,ndstep,ndata, ierr
  real dum
  character(len=256)  :: restart_dir, restart_fname
  character(len=1024) :: dest
  character(len=32)   :: diag_step_start_str
  character(len=32)   :: history_filename, data1d_filename

!! write out two copies of history.out and data1d.out for restart from a crash
     if(mype==0 .and. istep<=mstep)then
        if(mod(irest,2)==0)then
           restart_dir="restart_dir1/"
        else
           restart_dir="restart_dir2/"
        endif

#ifdef ADIOS2

        ! Copy the history.bp file
        dest = trim(adjustl(restart_dir))//"history_restart.bp"//char(0)
        call adios2_restart_hist (dest)
        write(gtcout,*)'write to ',trim(dest)
        
        ! Copy the data1d.bp file
        dest = trim(adjustl(restart_dir))//"data1d_restart.bp"//char(0)
        call adios2_restart_data1d (dest)
        write(gtcout,*)'write to ',trim(dest)

! #else
#endif

        ! history.out
        open(333,file=trim(restart_dir)//'history_restart.out',status='replace')
        rewind(iodiag)
        read(iodiag,101)mstepfinal
        
        ! total no. of diagnosis steps expected in the file == mstepfinal?
        ! mstep = no. of steps in this session
        ! ndiag = frequency of diagnosis 
        ! istep = current step index (starting from 1)
        ! ndstep (local) = no. of steps currently in file
        !   (previous steps + steps written in this session)
        ndstep=mstepfinal-mstep/ndiag+istep/ndiag
        write(333,101)ndstep
        do i=1,5
           read(iodiag,101)ii(i)
           write(333,101)ii(i)
        enddo

        ! this reads ndsteps of all the 3 arrays: partdata, fieldtime, fieldmode(1), fieldmode(2)
        ndata=ii(1)*ii(2)+ii(3)*(2*ii(4)+ii(5))
        do i=0,ndata*ndstep
           read(iodiag,102)dum
           write(333,102)dum
        enddo
        close(333)
        write(gtcout,*)'write to ',trim(restart_dir)//'history_restart.out'

        ! data1d.out
        open(444,file=trim(restart_dir)//'data1d_restart.out',status='replace')
        rewind(iodata1d)
        read(iodata1d,101)mstepfinal
        ndstep=mstepfinal-mstep/ndiag+istep/ndiag
        write(444,101)ndstep
        do i=1,6
           read(iodata1d,101)ii(i)
           write(444,101)ii(i)
        enddo


        ndata=ii(1)*(ii(2)*ii(4)+ii(5)*ii(6))
        do i=1,ndata*ndstep
           read(iodata1d,102)dum
           write(444,102)dum
        enddo
        close(444)
        write(gtcout,*)'write to ',trim(restart_dir)//'data1d_restart.out'

! end ifdef ADIOS2
! #endif

! save current location of restart files
        open(555,file="FileExit.out",status="replace")
        write(555,"(A9,i1)")"FileExit=",FileExit
        write(555,"(A9,i5)")"irest   =",irest+1
        if(mod(irest,2)==0)then
           write(555,"(A12)")"restart_dir1"
        else
           write(555,"(A12)")"restart_dir2"
        endif
        close(555)

        if(istep==mstep)then
! #ifndef ADIOS2
           close(iodiag)
           close(iodata1d)
! #endif
        endif

     endif

     irest=irest+1 !!XY rolling rstart
101  format(i6)
102  format(e13.6)

end subroutine restart_hist


#ifdef ADIOS2
  !> -------------------------------------------------------------------------
  !> Copy history.bp into history_restart.bp
  !> -------------------------------------------------------------------------
  subroutine adios2_restart_hist(dest)
    use adios2_common
    use adios2_history_diagnostics
    use adios2
    use mpi
    implicit none

    character(len=*), intent(in)  :: dest
    type(adios2_engine)           :: hist_engine, r_engine
    type(t_history_diagnosis)     :: t_hist_restart
    integer                       :: ierr, dummy_nsteps_in_file

    ! Open the history file for readin'
    call adios2_open (hist_engine, adios2_hist_restart_io, "history.bp", &
      adios2_mode_read, mpi_comm_self, ierr)

    ! Open the restart file for writin'
    call adios2_open (r_engine, adios2_hist_restart_io, trim(adjustl(dest)), &
      adios2_mode_write, mpi_comm_self, ierr)

    ! Call the copy subroutine
    ! We dont need to know the no. of steps in the history file. So,
    ! just pass a dummy intent(out) variable for that
    call copy_data (hist_engine, r_engine, adios2_hist_restart_io, &
      t_hist_restart, dummy_nsteps_in_file)

    ! Close the history and history_restart files
    call adios2_close (hist_engine, ierr)
    call adios2_close (r_engine,    ierr)

    ! Clear out the io object for reuse
    call adios2_remove_all_variables (adios2_hist_restart_io, ierr)

  end subroutine adios2_restart_hist


  !> -------------------------------------------------------------------------
  !> Copy data1d.bp into data1d_restart.bp
  !> -------------------------------------------------------------------------
  subroutine adios2_restart_data1d(dest)
    use adios2_common
    use adios2_data1d_diagnostics
    use adios2
    use mpi
    implicit none
    
    character(len=*), intent(in)  :: dest
    type(adios2_engine)           :: data1d_engine, r_engine
    type(t_data1d_diagnosis)      :: t_data1d_restart
    integer                       :: ierr

    ! Open the data1d file for readin'
    call adios2_open (data1d_engine, adios2_data1d_restart_io, "data1d.bp", &
      adios2_mode_read, mpi_comm_self, ierr)

    ! Open the restart file for writin'
    call adios2_open (r_engine, adios2_data1d_restart_io, trim(adjustl(dest)), &
      adios2_mode_write, mpi_comm_self, ierr)

    ! Call the copy subroutine
    ! We dont need to know the no. of steps in the data1d file. So,
    ! just pass a dummy intent(out) variable for that
    call copy_data (data1d_engine, r_engine, adios2_data1d_restart_io, t_data1d_restart)

    ! Close the data1d and data1d_restart files
    call adios2_close (data1d_engine, ierr)
    call adios2_close (r_engine,    ierr)

    ! Clear out the io object for reuse
    call adios2_remove_all_variables (adios2_data1d_restart_io, ierr)

  end subroutine adios2_restart_data1d
#endif

