#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_NOTOK_MDF(action) if (status/=0) then; TRACEBACK; action; call MDF_CLose(fid,status); status=1; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
#include "tm5.inc"
!----------------------------------------------------------------------------
!                  TM5                                                      !
!----------------------------------------------------------------------------
!BOP
!
! !MODULE:  Restart
!
! !DESCRIPTION: Write and read restart files.
!\\
!\\
! !INTERFACE: 
!
module Restart
  !
  ! !USES:
  !  
  use GO     , only : gol, goPr, goErr
  use dims   , only : nregions
  
  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public  ::  Restart_Init
  public  ::  Restart_Done
  public  ::  Restart_Save
  public  ::  Restart_Write
  public  ::  Restart_Read
  public  ::  rs_write
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'Restart'
  character(len=256)           ::  rs_write_dir
  logical                      ::  rs_write
  logical                      ::  rs_write_extra
  integer                      ::  rs_write_extra_dhour, rs_write_extra_hour
  integer                      ::  fid ! file id for IF_NOTOK_MDF macro
  !
  ! !REVISION HISTORY: 
  !    8 Apr 2011 - P. Le Sager - Close MDF file if error occurs. This is
  !                 needed for mpi_abort not to hang. See TM5_MPI_Abort in
  !                 partools, and remarks below. Made IF_NOTOK_MDF macro for
  !                 that purpose.
  !   28 Apr 2011 - P. Le Sager - Read method : handle restart file with extra
  !                 tracers. 
  !
  ! !REMARKS:
  ! (1) when an error occurs when accessing MDF files, you should first close
  !   the file before returning. The IF_NOTOK_MDF macro takes care of that.
  !   The only thing you need is to call it like that :
  !   
  !      IF_NOTOK_MDF(fid=xxxx)
  !
  !   where you replace xxxx with the integer id (file handler) of the file
  !   you are accessing. 
  ! 
  !EOP
  !------------------------------------------------------------------------

contains


  ! ================================================================
  
  
  subroutine Restart_Init( status )
  
    use GO         , only : TrcFile, Init, Done, ReadRc
    use global_data, only : rcfile
    use global_data, only : outdir
    use meteodata  , only : lli
  
    ! --- in/out -------------------------------
    
    integer, intent(out)    ::  status
    
    ! --- const --------------------------------
    
    character(len=*), parameter  ::  rname = 'Restart_Init'
    
    ! --- local --------------------------------
    
    type(TrcFile)      ::  rcF
    
    ! --- begin --------------------------------

    !
    ! read settings from rcfile
    !

    ! open:
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    ! write restart files at all ?
    call ReadRc( rcF, 'restart.write', rs_write, status, default=.false. )
    IF_ERROR_RETURN(status=1)
    
    ! further settings ...
    if ( rs_write ) then

      ! output directory:
      call ReadRc( rcF, 'restart.write.dir', rs_write_dir, status, default=outdir )
      IF_ERROR_RETURN(status=1)

      ! extra restart files ?
      call ReadRc( rcF, 'restart.write.extra', rs_write_extra, status, default=.false. )
      IF_ERROR_RETURN(status=1)
      if ( rs_write_extra ) then
        call ReadRc( rcF, 'restart.write.extra.hour', rs_write_extra_hour, status, default=0 )
        IF_ERROR_RETURN(status=1)
        call ReadRc( rcF, 'restart.write.extra.dhour', rs_write_extra_dhour, status, default=24 )
        IF_ERROR_RETURN(status=1)
      end if
      
    end if  ! write restart files

    ! close:
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)
    
    !
    ! done
    !
    
    ! ok
    status = 0
  
  end subroutine Restart_Init
  
  
  ! ***
  

  subroutine Restart_Done( status )
  
    ! --- in/out -------------------------------
    
    integer, intent(out)    ::  status
    
    ! --- const --------------------------------
    
    character(len=*), parameter  ::  rname = 'Restart_Done'
    
    ! --- begin --------------------------------
    
    ! nothing to be done ...
    
    ! ok
    status = 0
  
  end subroutine Restart_Done
  
  
  ! ***
  

  subroutine Restart_Save( status, extra, isfirst )
  
    use dims, only : idate
  
    ! --- in/out -------------------------------
    
    integer, intent(out)            ::  status
    logical, intent(in), optional   ::  extra
    logical, intent(in), optional   ::  isfirst
    
    ! --- const --------------------------------
    
    character(len=*), parameter  ::  rname = 'Restart_Save'
    
    ! --- local --------------------------------
    
    logical  ::  is_extra
    
    ! --- begin --------------------------------
    
    ! options ...
    is_extra = .false.
    if ( present(extra) ) is_extra = extra
    
    ! write restart files at all ?
    if ( rs_write ) then
    
      ! end or extra ?
      if ( is_extra ) then

        ! save extra restart files ?
        if ( rs_write_extra ) then

          ! every hour+n*dhour only :
          if (  modulo( idate(4) - rs_write_extra_hour, rs_write_extra_dhour ) == 0 .and. &
                all( idate(5:6) == 0 ) ) then

            ! write restart file for this time:
            call Restart_Write( status, isfirst=isfirst )
            IF_NOTOK_RETURN(status=1)

          end if  ! for this hour

        end if   ! extra restart files ?

      else

        ! write restart file :
        call Restart_Write( status, isfirst=isfirst )
        IF_NOTOK_RETURN(status=1)

      end if  ! not extra
    
    end if  ! write at all

    ! ok
    status = 0
  
  end subroutine Restart_Save
  
  
  ! ***
  
  subroutine Restart_FileName( region, fname, status, key, dir, isfirst )
  
    use dims       , only : idate
    use global_data, only : outdir
    use meteodata  , only : lli
  
    ! --- in/out -------------------------------
    
    integer, intent(in)             ::  region
    character(len=*), intent(out)   ::  fname
    integer, intent(out)            ::  status
 
    character(len=*), intent(in), optional  ::  dir    
    character(len=*), intent(in), optional  ::  key
    logical, intent(in), optional           ::  isfirst
    
    ! --- local -------------------------------
    
    character(len=*), parameter  ::  rname = 'Restart_FileName'
    character(len=256)           ::  adir
    character(len=32)            ::  akey
    
    ! --- begin --------------------------------
 
    ! destination directory:
    adir = trim(outdir)
    if ( present(dir) ) adir = trim(dir)

    ! extra key, for example '_x' to denote that 
    ! a restart file was dumped after process 'x':
    akey = ''
    if ( present(key) )  akey = trim(key)
    
    ! if this is the initial time, add an extra key to avoid
    ! that the restart file for this hour from the previous
    ! run is overwritten:
    if ( present(isfirst) ) then
      if ( isfirst ) akey = trim(akey)//'_initial'
    end if

    ! write filename:
    write (fname,'(a,"/TM5_restart_",i4.4,2i2.2,"_",2i2.2,"_",a,a,".nc")') &
                trim(adir), idate(1:5), trim(lli(region)%name), trim(akey)

    ! ok
    status = 0
  
  end subroutine Restart_FileName
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  Restart_Write
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Restart_Write( status, key, region, isfirst )

    use GO          , only : Get
    use dims        , only : nregions, at, bt
!    use dims        , only : nlon360, nlat180, iglbsfc
    use chem_param  , only : ntracet, ntrace_chem, ntrace, names
    use partools    , only : localComm
    use partools    , only : myid, root
    use partools    , only : ntracetloc, offsetn
    use partools    , only : lmloc, offsetl
    use partools    , only : MPI_INFO_NULL
    use partools    , only : previous_par
    use tracer_data , only : mass_dat, chem_dat
#ifdef with_tendencies
    use tm5_tendency, only : plc_ntr, plc_trname
    use tm5_tendency, only : plc_npr, plc_prname
    use tracer_data , only : plc_dat
#endif
    use meteodata   , only : lli, levi
    use meteodata   , only : sp_dat, phlb_dat, m_dat
    !use meteodata   , only : sp1_dat, sp2_dat                 ! debug
    !use meteodata   , only : mfu_dat, mfv_dat, mfw_dat        ! debug
    !use meteodata   , only : pu_dat, pv_dat, pw_dat           ! debug
    !use meteodata   , only : sshf_dat, slhf_dat
    use MDF         , only : MDF_Create, MDF_EndDef, MDF_Close
    use MDF         , only : MDF_Var_Par_Access, MDF_INDEPENDENT
    use MDF         , only : MDF_Def_Dim, MDF_Def_Var
    use MDF         , only : MDF_Put_Att
    use MDF         , only : MDF_Put_Var
    use MDF         , only : MDF_REPLACE, MDF_NETCDF, MDF_NETCDF4
    use MDF         , only : MDF_FLOAT, MDF_DOUBLE, MDF_CHAR
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,          intent(out)           ::  status
    !
    ! !INPUT PARAMETERS:
    !    
    character(len=*), intent(in), optional  ::  key
    integer,          intent(in), optional  ::  region
    logical,          intent(in), optional  ::  isfirst
    !
    ! !REVISION HISTORY: 
    !    8 Apr 2011 - P. Le Sager - use IF_NOTOK_MDF macro
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
        
    ! --- local --------------------------------
    
    character(len=*), parameter  ::  rname = 'Restart_Write'
    
    integer               ::  imr, jmr, lmr
    integer               ::  iregion
    character(len=256)    ::  fname
    integer               ::  ftype
    integer               ::  ncid
    integer               ::  dimid_lon, dimid_lat, dimid_lev, dimid_hlev
    integer               ::  dimid_lon_sfc, dimid_lat_sfc
    integer               ::  dimid_trace, dimid_trace_transp, dimid_trace_chem
    integer               ::  dimid_name
    integer               ::  varid
    integer               ::  varid_sp, varid_ph, varid_m, varid_at, varid_bt
    !integer               ::  varid_sp1, varid_sp2               ! debug
    !integer               ::  varid_mfu, varid_mfv, varid_mfw    ! debug
    !integer               ::  varid_pu, varid_pv, varid_pw       ! debug
    !integer               ::  varid_sshf, varid_slhf
    integer               ::  varid_names
    integer               ::  varid_rm
#ifdef slopes
    integer               ::  varid_rxm, varid_rym, varid_rzm
#endif
    integer               ::  varid_rmc
#ifdef with_tendencies
    integer               ::  varid_plc(plc_ntr,plc_npr)
    integer               ::  itr, ipr
    integer               ::  time6(6)
#endif
    integer               ::  rtype
    
    ! --- begin --------------------------------
    
    write (gol,'("write restart file(s) ...")'); call goPr
    
    ! loop over regions:
    do iregion = 1, nregions
    
      ! only selected region ?
      if ( present(region) ) then
        if ( iregion /= region ) cycle
      end if
    
      ! grid size
      imr = lli(iregion)%nlon
      jmr = lli(iregion)%nlat
      lmr = levi%nlev
    
      ! name of restart file
      call Restart_FileName( iregion, fname, status, key=key, dir=rs_write_dir, isfirst=isfirst )
      IF_NOTOK_RETURN(status=1)
      
      write (gol,'("  destination : ",a)') trim(fname); call goPr

      ! o open netcdf file
      
#ifdef MPI      
      ! overwrite existing files (clobber), provide MPI stuff:
      call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, ncid, status, &
                          mpi_comm=localComm, mpi_info=MPI_INFO_NULL )
      if (status/=0) then
        write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr
        write (gol,'("MDF module not compiled with netcdf4 support ?")'); call goErr
        TRACEBACK; status=1; return
      end if
#else
      ! overwrite existing files (clobber)
      call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, ncid, status )
      IF_NOTOK_RETURN(status=1)
#endif
      
      ! o define dimensions

      call MDF_Def_Dim( ncid, 'lon', lli(iregion)%nlon, dimid_lon, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Def_Dim( ncid, 'lat', lli(iregion)%nlat, dimid_lat, status )
      IF_NOTOK_MDF(fid=ncid)
! These Dims are not used, so commented (pls, 30-03-2011)
!      call MDF_Def_Dim( ncid, 'lon_sfc', lli(iglbsfc)%nlon, dimid_lon_sfc, status )
!      IF_NOTOK_MDF(fid=ncid)
!      call MDF_Def_Dim( ncid, 'lat_sfc', lli(iglbsfc)%nlat, dimid_lat_sfc, status )
!      IF_NOTOK_MDF(fid=ncid)
      call MDF_Def_Dim( ncid, 'lev', levi%nlev, dimid_lev, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Def_Dim( ncid, 'hlev', levi%nlev+1, dimid_hlev, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Def_Dim( ncid, 'trace_transp', ntracet, dimid_trace_transp, status )
      IF_NOTOK_MDF(fid=ncid)
      if ( ntrace_chem > 0 ) then
        call MDF_Def_Dim( ncid, 'trace_chem', ntrace_chem, dimid_trace_chem, status )
        IF_NOTOK_MDF(fid=ncid)
      else
        dimid_trace_chem = -1
      end if
      call MDF_Def_Dim( ncid, 'trace', ntrace, dimid_trace, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Def_Dim( ncid, 'name', len(names(1)), dimid_name, status )
      IF_NOTOK_MDF(fid=ncid)

      ! o define variables
      
      select case ( kind(m_dat(iregion)%data) )
        case ( 4 ) ; rtype = MDF_FLOAT
        case ( 8 ) ; rtype = MDF_DOUBLE
        case default
          write (gol,'("unsupported real kind : ",i6)') kind(m_dat(iregion)%data)
          TRACEBACK; status=1; return
      end select

      ! surface pressure
      call MDF_Def_Var( ncid, 'sp', rtype, &
                               (/dimid_lon,dimid_lat/), varid, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'long_name', 'surface pressure', status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'unit', 'Pa', status )
      IF_NOTOK_MDF(fid=ncid)
      varid_sp = varid

      ! at, bt coefficients for hybrid grid
      call MDF_Def_Var( ncid, 'at', rtype, (/dimid_hlev/), varid, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'long_name', 'hybrid grid a_t coefficient', status )
      IF_NOTOK_MDF(fid=ncid)
      varid_at = varid

      call MDF_Def_Var( ncid, 'bt', rtype, (/dimid_hlev/), varid, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'long_name', 'hybrid grid b_t coefficient', status )
      IF_NOTOK_MDF(fid=ncid)
      varid_bt = varid

      !! debug ...
      !MDF_Def_Var( ncid, 'sp1', rtype, (/dimid_lon,dimid_lat/), varid_sp1, status )
      !IF_NOTOK_MDF(fid=ncid)
      !MDF_Def_Var( ncid, 'sp2', rtype, (/dimid_lon,dimid_lat/), varid_sp2, status )
      !IF_NOTOK_MDF(fid=ncid)

      ! half level pressure
      call MDF_Def_Var( ncid, 'ph', rtype, &
                               (/dimid_lon,dimid_lat,dimid_hlev/), varid, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'long_name', 'half level pressure', status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'unit', 'Pa', status )
      IF_NOTOK_MDF(fid=ncid)
      varid_ph = varid
      
      ! air mass
      call MDF_Def_Var( ncid, 'm', rtype, &
                               (/dimid_lon,dimid_lat,dimid_lev/), varid, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'long_name', 'air mass', status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'unit', 'kg', status )
      IF_NOTOK_MDF(fid=ncid)
      varid_m = varid
      
      !! debug
      !call MDF_Def_Var( ncid, 'mfu', rtype, (/dimid_lon,dimid_lat,dimid_lev/), varid_mfu, status )
      !IF_NOTOK_MDF(fid=ncid)
      !call MDF_Def_Var( ncid, 'mfv', rtype, (/dimid_lon,dimid_lat,dimid_lev/), varid_mfv, status )
      !IF_NOTOK_MDF(fid=ncid)
      !call MDF_Def_Var( ncid, 'mfw', rtype, (/dimid_lon,dimid_lat,dimid_lev/), varid_mfw, status )
      !IF_NOTOK_MDF(fid=ncid)
      !call MDF_Def_Var( ncid, 'pu', rtype, (/dimid_lon,dimid_lat,dimid_lev/), varid_pu, status )
      !IF_NOTOK_MDF(fid=ncid)
      !call MDF_Def_Var( ncid, 'pv', rtype, (/dimid_lon,dimid_lat,dimid_lev/), varid_pv, status )
      !IF_NOTOK_MDF(fid=ncid)
      !call MDF_Def_Var( ncid, 'pw', rtype, (/dimid_lon,dimid_lat,dimid_lev/), varid_pw, status )
      !IF_NOTOK_MDF(fid=ncid)
      
      !! accumulated surface fluxes
      !!
      !call MDF_Def_Var( ncid, 'slhf', rtype, (/dimid_lon_sfc,dimid_lat_sfc/), varid, status )
      !IF_NOTOK_MDF(fid=ncid)
      !call MDF_Put_Att( ncid, varid, 'long_name', 'surface latent heat flux', status )
      !IF_NOTOK_MDF(fid=ncid)
      !call MDF_Put_Att( ncid, varid, 'unit', 'W/m2', status )
      !IF_NOTOK_MDF(fid=ncid)
      !varid_slhf = varid
      !!
      !call MDF_Def_Var( ncid, 'sshf', rtype, (/dimid_lon_sfc,dimid_lat_sfc/), varid, status )
      !IF_NOTOK_MDF(fid=ncid)
      !call MDF_Put_Att( ncid, varid, 'long_name', 'surface sensible heat flux', status )
      !IF_NOTOK_MDF(fid=ncid)
      !call MDF_Put_Att( ncid, varid, 'unit', 'W/m2', status )
      !IF_NOTOK_MDF(fid=ncid)
      !varid_sshf = varid
      
      ! tracer names
      call MDF_Def_Var( ncid, 'names', MDF_CHAR, (/dimid_name,dimid_trace/), varid, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'long_name', 'tracer names', status )
      IF_NOTOK_MDF(fid=ncid)
      varid_names = varid
      
      ! tracer mass
      call MDF_Def_Var( ncid, 'rm', rtype, &
                               (/dimid_lon,dimid_lat,dimid_lev,dimid_trace_transp/), varid, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'long_name', 'transported tracer mass', status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'unit', 'kg', status )
      IF_NOTOK_MDF(fid=ncid)
      varid_rm = varid
      
      ! tracer mass slopes:
#ifdef slopes        
      call MDF_Def_Var( ncid, 'rxm', rtype, &
                               (/dimid_lon,dimid_lat,dimid_lev,dimid_trace_transp/), varid, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'long_name', 'tracer mass slope in x direction', status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'unit', 'kg/(half cell)', status )
      IF_NOTOK_MDF(fid=ncid)
      varid_rxm = varid
      
      call MDF_Def_Var( ncid, 'rym', rtype, &
                               (/dimid_lon,dimid_lat,dimid_lev,dimid_trace_transp/), varid, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'long_name', 'tracer mass slope in y direction', status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'unit', 'kg/(half cell)', status )
      IF_NOTOK_MDF(fid=ncid)
      varid_rym = varid
      
      call MDF_Def_Var( ncid, 'rzm', rtype, &
                               (/dimid_lon,dimid_lat,dimid_lev,dimid_trace_transp/), varid, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'long_name', 'tracer mass slope in z direction', status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Put_Att( ncid, varid, 'unit', 'kg/(half cell)', status )
      IF_NOTOK_MDF(fid=ncid)
      varid_rzm = varid
#endif

      ! non-transported tracers:
      if ( ntrace_chem > 0 ) then 
        call MDF_Def_Var( ncid, 'rmc', rtype, &
                                 (/dimid_lon,dimid_lat,dimid_lev,dimid_trace_chem/), varid, status )
        IF_NOTOK_MDF(fid=ncid)
        call MDF_Put_Att( ncid, varid, 'long_name', 'non-transported tracer mass', status )
        IF_NOTOK_MDF(fid=ncid)
        call MDF_Put_Att( ncid, varid, 'unit', 'kg', status )
        IF_NOTOK_MDF(fid=ncid)
        varid_rmc = varid
      end if
      
#ifdef with_tendencies
      ! production, loss, and concentration:
      do itr = 1, plc_ntr
        do ipr = 1, plc_npr
          ! define netcdf variable:
          call MDF_Def_Var( ncid, trim(plc_trname(itr))//'_'//trim(plc_prname(ipr)), rtype, &
                                   (/dimid_lon,dimid_lat,dimid_lev/), varid, status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Put_Att( ncid, varid, 'long_name', 'chemical tendency', status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Put_Att( ncid, varid, 'unit', trim(plc_dat(region,itr,ipr)%unit), status )
          IF_NOTOK_MDF(fid=ncid)
          ! extract time as 6 integers:
          call Get( plc_dat(region,itr,ipr)%t, time6=time6 )
          ! add time attribute:
          call MDF_Put_Att( ncid, varid, 'time', time6, status )
          IF_NOTOK_MDF(fid=ncid)
          ! store variable id:
          varid_plc(itr,ipr) = varid
        end do
      end do
#endif

      !  o end definition mode

      call MDF_EndDef( ncid, status )
      IF_NOTOK_MDF(fid=ncid)
      
      ! o write variables

      ! set independent data mode (not all processes are writing):
      call MDF_Var_Par_Access( ncid, varid_sp  , MDF_INDEPENDENT, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Var_Par_Access( ncid, varid_ph  , MDF_INDEPENDENT, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Var_Par_Access( ncid, varid_m   , MDF_INDEPENDENT, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Var_Par_Access( ncid, varid_at   , MDF_INDEPENDENT, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Var_Par_Access( ncid, varid_bt   , MDF_INDEPENDENT, status )
      IF_NOTOK_MDF(fid=ncid)
      !call MDF_Var_Par_Access( ncid, varid_slhf, MDF_INDEPENDENT, status )
      !IF_NOTOK_MDF(fid=ncid)
      !call MDF_Var_Par_Access( ncid, varid_sshf, MDF_INDEPENDENT, status )
      !IF_NOTOK_MDF(fid=ncid)

      ! write meteo etc from root:
      if ( myid == root ) then

        ! surface pressure
        call MDF_Put_Var( ncid, varid_sp, sp_dat(iregion)%data(1:imr,1:jmr,1), status )
        IF_NOTOK_MDF(fid=ncid)

        !! debug ...
        !call MDF_Put_Var( ncid, varid_sp1, sp1_dat(iregion)%data(1:imr,1:jmr,1), status )
        !IF_NOTOK_MDF(fid=ncid)
        !call MDF_Put_Var( ncid, varid_sp2, sp2_dat(iregion)%data(1:imr,1:jmr,1), status )
        !IF_NOTOK_MDF(fid=ncid)
        
        ! half level pressure
        call MDF_Put_Var( ncid, varid_ph, phlb_dat(iregion)%data(1:imr,1:jmr,1:lmr+1), status )
        IF_NOTOK_MDF(fid=ncid)
        
        ! at, bt coefficients
        call MDF_Put_Var( ncid, varid_at, at(1:lmr+1), status )
        IF_NOTOK_MDF(fid=ncid)
        
        call MDF_Put_Var( ncid, varid_bt, bt(1:lmr+1), status )
        IF_NOTOK_MDF(fid=ncid)
        
        ! air mass
        call MDF_Put_Var( ncid, varid_m, m_dat(iregion)%data(1:imr,1:jmr,1:lmr), status )
        IF_NOTOK_MDF(fid=ncid)
        
        !! debug ...
        !call MDF_Put_Var( ncid, varid_mfu, mfu_dat(iregion)%data(1:imr,1:jmr,1:lmr), status )
        !IF_NOTOK_MDF(fid=ncid)
        !call MDF_Put_Var( ncid, varid_mfv, mfv_dat(iregion)%data(1:imr,1:jmr,1:lmr), status )
        !IF_NOTOK_MDF(fid=ncid)
        !call MDF_Put_Var( ncid, varid_mfw, mfw_dat(iregion)%data(1:imr,1:jmr,1:lmr), status )
        !IF_NOTOK_MDF(fid=ncid)
        !call MDF_Put_Var( ncid, varid_pu, pu_dat(iregion)%data(1:imr,1:jmr,1:lmr), status )
        !IF_NOTOK_MDF(fid=ncid)
        !call MDF_Put_Var( ncid, varid_pv, pv_dat(iregion)%data(1:imr,1:jmr,1:lmr), status )
        !IF_NOTOK_MDF(fid=ncid)
        !call MDF_Put_Var( ncid, varid_pw, pw_dat(iregion)%data(1:imr,1:jmr,1:lmr), status )
        !IF_NOTOK_MDF(fid=ncid)

        !! surface latent heat flux; global surface field !
        !call MDF_Put_Var( ncid, varid_slhf, slhf_dat(iglbsfc)%data(1:nlon360,1:nlat180,1), status )
        !IF_NOTOK_MDF(fid=ncid)
        !
        !! surface sensible heat flux; global surface field !
        !call MDF_Put_Var( ncid, varid_sshf, sshf_dat(iglbsfc)%data(1:nlon360,1:nlat180,1), status )
        !IF_NOTOK_MDF(fid=ncid)
        
        ! tracer names:
        call MDF_Put_Var( ncid, varid_names, names, status )
        IF_NOTOK_MDF(fid=ncid)

      end if  ! root
      
      ! set independent data mode (not all processes are writing):
      call MDF_Var_Par_Access( ncid, varid_rm , MDF_INDEPENDENT, status )
      IF_NOTOK_MDF(fid=ncid)
#ifdef slopes   
      call MDF_Var_Par_Access( ncid, varid_rxm, MDF_INDEPENDENT, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Var_Par_Access( ncid, varid_rym, MDF_INDEPENDENT, status )
      IF_NOTOK_MDF(fid=ncid)
      call MDF_Var_Par_Access( ncid, varid_rzm, MDF_INDEPENDENT, status )
      IF_NOTOK_MDF(fid=ncid)
#endif

      ! write transported tracers, distributed over tracers or levels:
      select case ( previous_par(iregion) )
      
        case ( 'tracer' )
        
          ! distributed over tracers; might be empty on this pe
          if ( ntracetloc > 0 ) then

            call MDF_Put_Var( ncid, varid_rm, &
                          mass_dat(iregion)%rm_t(1:imr,1:jmr,1:lmr,1:ntracetloc), status, &
                          start=(/1,1,1,offsetn+1/), &
                          count=(/imr,jmr,lmr,ntracetloc/) )
            IF_NOTOK_MDF(fid=ncid)

#ifdef slopes        
            call MDF_Put_Var( ncid, varid_rxm, &
                            mass_dat(iregion)%rxm_t(1:imr,1:jmr,1:lmr,1:ntracetloc), status, &
                          start=(/1,1,1,offsetn+1/), &
                          count=(/imr,jmr,lmr,ntracetloc/) )
            IF_NOTOK_MDF(fid=ncid)

            call MDF_Put_Var( ncid, varid_rym, &
                            mass_dat(iregion)%rym_t(1:imr,1:jmr,1:lmr,1:ntracetloc), status, &
                          start=(/1,1,1,offsetn+1/), &
                          count=(/imr,jmr,lmr,ntracetloc/) )
            IF_NOTOK_MDF(fid=ncid)

            call MDF_Put_Var( ncid, varid_rzm, &
                            mass_dat(iregion)%rzm_t(1:imr,1:jmr,1:lmr,1:ntracetloc), status, &
                          start=(/1,1,1,offsetn+1/), &
                          count=(/imr,jmr,lmr,ntracetloc/) )
            IF_NOTOK_MDF(fid=ncid)
#endif
          end if   ! ntraceloc > 0
          
        case ( 'levels' )
        
          ! distributed over levels; might be empty on this pe
          if ( lmloc > 0 ) then

            call MDF_Put_Var( ncid, varid_rm, &
                            mass_dat(iregion)%rm_k(1:imr,1:jmr,1:lmloc,1:ntracet), status, &
                            start=(/1,1,offsetl+1,1/), &
                            count=(/imr,jmr,lmloc,ntracet/) )
            IF_NOTOK_MDF(fid=ncid)

#ifdef slopes        
            call MDF_Put_Var( ncid, varid_rxm, &
                            mass_dat(iregion)%rxm_k(1:imr,1:jmr,1:lmloc,1:ntracet), status, &
                            start=(/1,1,offsetl+1,1/), &
                            count=(/imr,jmr,lmloc,ntracet/) )
            IF_NOTOK_MDF(fid=ncid)

            call MDF_Put_Var( ncid, varid_rym, &
                            mass_dat(iregion)%rym_k(1:imr,1:jmr,1:lmloc,1:ntracet), status, &
                            start=(/1,1,offsetl+1,1/), &
                            count=(/imr,jmr,lmloc,ntracet/) )
            IF_NOTOK_MDF(fid=ncid)

            call MDF_Put_Var( ncid, varid_rzm, &
                            mass_dat(iregion)%rzm_k(1:imr,1:jmr,1:lmloc,1:ntracet), status, &
                            start=(/1,1,offsetl+1,1/), &
                            count=(/imr,jmr,lmloc,ntracet/) )
            IF_NOTOK_MDF(fid=ncid)
#endif
          end if   ! ntraceloc > 0
          
        case default

          write (gol,'("unsupported distribution:")'); call goPr
          write (gol,'("  region         : ",i2)') iregion ;call goPr
          write (gol,'("  previous_par   : ",a)') previous_par(iregion) ;call goPr
          TRACEBACK; status=1; return

      end select
      
      ! write non-transported tracers levels on this pe:
      if ( (ntrace_chem > 0) .and. (lmloc > 0) ) then

        ! set independent data mode (not all processes are writing):
        call MDF_Var_Par_Access( ncid, varid_rmc, MDF_INDEPENDENT, status )
        IF_NOTOK_MDF(fid=ncid)

      
        call MDF_Put_Var( ncid, varid_rmc, &
                        chem_dat(iregion)%rm_k(1:imr,1:jmr,1:lmloc,ntracet+1:ntracet+ntrace_chem), status, &
                        start=(/1,1,offsetl+1,1/), &
                        count=(/imr,jmr,lmloc,ntrace_chem/) )
        IF_NOTOK_MDF(fid=ncid)
        
      end if   ! lmloc > 0
      
#ifdef with_tendencies
      ! set independent data mode (not all processes are writing):
      do itr = 1, plc_ntr
        do ipr = 1, plc_npr
          call MDF_Var_Par_Access( ncid, varid_plc(itr,ipr) , MDF_INDEPENDENT, status )
          IF_NOTOK_MDF(fid=ncid)
        end do
      end do

      ! write production/loss/concentration levels on this pe:
      if ( lmloc > 0 ) then
        
        do itr = 1, plc_ntr
          do ipr = 1, plc_npr
            call MDF_Put_Var( ncid, varid_plc(itr,ipr), &
                            plc_dat(region,itr,ipr)%rm_k(1:imr,1:jmr,1:lmloc), status, &
                            start=(/1,1,offsetl+1/), count=(/imr,jmr,lmloc/) )
            IF_NOTOK_MDF(fid=ncid)
          end do
        end do
        
      end if
#endif

      ! o close file
      call MDF_Close( ncid, status )
      IF_NOTOK_RETURN(status=1)

    end do   ! regions

    !write (gol,'("  ok")'); call goPr
    
    ! ok
    status = 0
    
  end subroutine Restart_Write
  !EOC
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  Restart_Read
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Restart_Read( status, &
                             region, &
                             surface_pressure, pressure, air_mass, surface_fluxes, &
                             tracer_mass, tendencies, megan_history, nox_pulsing, &
                             override_istart)
    !
    ! !USES:
    !
    use GO          , only : TrcFile, Init, Done, ReadRc
    use GO          , only : Set
    use GO          , only : goMatchValue
    use dims        , only : nregions
    use dims        , only : istart, idate, idatei
!    use dims        , only : nlon360, nlat180, iglbsfc
    use chem_param  , only : ntracet, ntrace_chem, ntrace
    use chem_param  , only : names, tracer_name_len
    use global_data , only : rcfile
    use partools    , only : localComm
    use partools    , only : myid, root
    use partools    , only : ntracetloc, offsetn, tracer_active, tracer_loc
    use partools    , only : lmloc, offsetl
    use partools    , only : MPI_INFO_NULL
    use partools    , only : previous_par
    use tracer_data , only : mass_dat, chem_dat
#ifdef with_tendencies
    use tm5_tendency, only : plc_ntr, plc_trname
    use tm5_tendency, only : plc_npr, plc_prname
    use tracer_data , only : plc_dat
#endif
    use meteodata   , only : lli, levi
    use meteodata   , only : sp_dat, phlb_dat, m_dat
    !use meteodata   , only : slhf_dat, sshf_dat
    use MDF         , only : MDF_Open,           MDF_Close,            MDF_Inquire_Dimension
    use MDF         , only : MDF_Inq_VarID,      MDF_Inquire_Variable, MDF_Inq_DimID
    use MDf         , only : MDF_Var_Par_Access, MDF_INDEPENDENT
    use MDF         , only : MDF_Get_Att
    use MDF         , only : MDF_Get_Var
    use MDF         , only : MDF_READ, MDF_NETCDF4
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)             ::  status
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in), optional    ::  region
    logical, intent(in), optional    ::  surface_pressure, pressure, air_mass, surface_fluxes
    logical, intent(in), optional    ::  tracer_mass, tendencies, megan_history, nox_pulsing
    logical, intent(in), optional    ::  override_istart
    !
    ! !REVISION HISTORY: 
    !    8 Apr 2011 - P. Le Sager - use IF_NOTOK_MDF macro
    !   28 Apr 2011 - P. Le Sager - Check on tracer availability in restart file.
    !                             - Allows for more tracers in restart file than needed
    !   10 May 2011 - P. Le Sager - Added deallocate statement to work with zoom regions
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
    ! --- const --------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/Restart_Read'
    
    ! --- local --------------------------------
    
    type(TrcFile)                   ::  rcF    
    character(len=256)              ::  rs_read_dir
    logical                         ::  exist
    logical                         ::  do_sp, do_ph, do_m, do_sflux, do_rm, do_plc, do_megan, do_pulse
    integer                         ::  imr, jmr, lmr
    integer                         ::  n, region1, region2
    character(len=256)              ::  fname
    integer                         ::  ncid
    integer                         ::  varid_sp
    integer                         ::  varid_ph
    integer                         ::  varid_m
    !integer                         ::  varid_slhf, varid_sshf
    integer                         ::  varid_names
    integer                         ::  iname
    character(len=tracer_name_len), allocatable ::  values_names(:)
    integer                         ::  itr, itr_loc, itr_file
    integer                         ::  varid_rm
    integer                         ::  ntracet_restart, dimid
    integer                         ::  shp(2)
#ifdef slopes
    integer                         ::  varid_rxm, varid_rym, varid_rzm
#endif
    integer                         ::  varid_rmc
#ifdef with_tendencies
    integer                         ::  varid_plc(plc_ntr,plc_npr)
    integer                         ::  itr, ipr
    integer                         ::  time6(6)
#endif

    
    ! --- begin --------------------------------
    
    write (gol,'("read restart file ...")'); call goPr
    
    ! correct istart ?
    if ( istart /= 33 .and. .not. present(override_istart) ) then
      write (gol,'("  skip; istart not 33 but ",i3)') istart; call goPr
      status=0; return
    end if
    
    ! init time ?
    if ( any( idate /= idatei ) ) then
      write (gol,'("  skip; idate not idatei but ",i4,5i2.2)') idate; call goPr
      status=0; return
    end if

    ! input directory:
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'restart.read.dir', rs_read_dir, status )
    IF_NOTOK_RETURN(status=1)
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

    ! region range:
    if ( present(region) ) then
      region1 = region
      region2 = region
    else
      region1 = 1
      region2 = nregions
    end if
    
    ! data sets:
    do_sp      = .false.  ;  if ( present(surface_pressure ) ) do_sp      = surface_pressure
    do_ph      = .false.  ;  if ( present(pressure         ) ) do_ph      = pressure
    do_m       = .false.  ;  if ( present(air_mass         ) ) do_m       = air_mass
    do_sflux   = .false.  ;  if ( present(surface_fluxes   ) ) do_sflux   = surface_fluxes
    do_rm      = .false.  ;  if ( present(tracer_mass      ) ) do_rm      = tracer_mass
    do_plc     = .false.  ;  if ( present(tendencies       ) ) do_plc     = tendencies
    do_megan   = .false.  ;  if ( present(megan_history    ) ) do_megan   = megan_history
    do_pulse   = .false.  ;  if ( present(nox_pulsing      ) ) do_pulse   = nox_pulsing
    
    ! sorry ..
    if ( do_sflux ) then
      write (gol,'("no surface fluxes in restart files until somebody")'); call goErr
      write (gol,'("has a good idea on what should be storred:")'); call goErr
      write (gol,'(" o global surface field (1x1 ?)")'); call goErr
      write (gol,'(" o zoom regions")'); call goErr
      write (gol,'(" o both")'); call goErr
      TRACEBACK; status=1; return
    end if
    
    ! loop over ns
    do n = region1, region2

      ! grid size
      imr = lli(n)%nlon
      jmr = lli(n)%nlat
      lmr = levi%nlev
    
      ! name of restart file
      call Restart_FileName( n, fname, status, dir=trim(rs_read_dir) )
      IF_NOTOK_RETURN(status=1)
      
      write (gol,'("  restore from ",a)') trim(fname); call goPr
    
      ! test ...
      inquire( file=fname, exist=exist )
      if ( .not. exist ) then
        write (gol,'("restart file not found : ",a)') trim(fname); call goErr
        TRACEBACK; status=1; return
      end if

      ! o open netcdf file

#if defined MPI && defined with_netcdf4_par      
      ! open existing file:
      call MDF_Open( trim(fname), MDF_NETCDF4, MDF_READ, ncid, status, &
                        mpi_comm=localComm, mpi_info=MPI_INFO_NULL )
      IF_NOTOK_RETURN(status=1)
#else
      ! open existing file:
      call MDF_Open( trim(fname), MDF_NETCDF4, MDF_READ, ncid, status )
      IF_NOTOK_RETURN(status=1)
#endif
      
      ! o define variables
      
      ! surface pressure
      if ( do_sp ) then
        call MDF_Inq_VarID( ncid, 'sp', varid_sp, status )
        IF_NOTOK_MDF(fid=ncid)
      end if
      
      ! half level pressure
      if ( do_ph ) then
        call MDF_Inq_VarID( ncid, 'ph', varid_ph, status )
        IF_NOTOK_MDF(fid=ncid)
      end if
      
      ! air mass
      if ( do_m ) then
        call MDF_Inq_VarID( ncid, 'm', varid_m, status )
        IF_NOTOK_MDF(fid=ncid)
      end if

      !! surface fluxes; only once is ok
      !if ( do_sflux .and. (n==region1) ) then
      !  call MDF_Inq_VarID( ncid, 'slhf', varid_slhf, status )
      !  IF_NOTOK_MDF(fid=ncid)
      !  call MDF_Inq_VarID( ncid, 'sshf', varid_sshf, status )
      !  IF_NOTOK_MDF(fid=ncid)
      !end if

      ! tracer mass
      if ( do_rm ) then

        call MDF_Inq_VarID( ncid, 'names', varid_names, status )
        if ( status /= 0 ) then
          write (gol,'("could not find variable `names` in restart file;")'); call goErr
          write (gol,'("  using an old restart file to initialize the model ?")'); call goErr
          status=1
        end if
        IF_NOTOK_MDF(fid=ncid)
        
        ! get dimension of "names"
        call MDF_Inquire_Variable( ncid, varid_names, status, shp=shp )
        IF_NOTOK_MDF(fid=ncid)

        ! get number of transported tracer in restart file
        call MDF_Inq_DimID( ncid, 'trace_transp', dimid, status )
        IF_NOTOK_MDF(fid=ncid)

        call MDF_Inquire_Dimension( ncid, dimid, status, length=ntracet_restart )
        IF_NOTOK_MDF(fid=ncid)

        ! various id
        call MDF_Inq_VarID( ncid, 'rm', varid_rm, status )
        IF_NOTOK_MDF(fid=ncid)

#ifdef slopes        
        call MDF_Inq_VarID( ncid, 'rxm', varid_rxm, status )
        IF_NOTOK_MDF(fid=ncid)

        call MDF_Inq_VarID( ncid, 'rym', varid_rym, status )
        IF_NOTOK_MDF(fid=ncid)

        call MDF_Inq_VarID( ncid, 'rzm', varid_rzm, status )
        IF_NOTOK_MDF(fid=ncid)
#endif

        ! read non-transported tracers levels on this pe:
        if ( ntrace_chem > 0 ) then
           call MDF_Inq_VarID( ncid, 'rmc', varid_rmc, status )
           IF_NOTOK_MDF(fid=ncid)
        end if
      end if
      
#ifdef with_tendencies
      ! production/loss/concentration
      if ( do_plc ) then
        do itr = 1, plc_ntr
          do ipr = 1, plc_npr
            call MDF_Inq_VarID( ncid, trim(plc_trname(itr))//'_'//trim(plc_prname(ipr)), varid_plc(itr,ipr), status )
            IF_NOTOK_MDF(fid=ncid)
          end do
        end do
      end if
#endif
      
      ! o read variables
    
      ! read meteo on each processor:

      ! surface pressure
      if ( do_sp ) then
        write (gol,'("  restore surface pressure ...")'); call goPr
        call MDF_Var_Par_Access( ncid, varid_sp, MDF_INDEPENDENT, status )
        IF_NOTOK_MDF(fid=ncid)
        call MDF_Get_Var( ncid, varid_sp, sp_dat(n)%data(1:imr,1:jmr,1), status )
        IF_NOTOK_MDF(fid=ncid)
      end if

      ! half level pressure
      if ( do_ph ) then
        write (gol,'("  restore half level pressure ...")'); call goPr
        call MDF_Var_Par_Access( ncid, varid_ph, MDF_INDEPENDENT, status )
        IF_NOTOK_MDF(fid=ncid)
        call MDF_Get_Var( ncid, varid_ph, phlb_dat(n)%data(1:imr,1:jmr,1:lmr+1), status )
        IF_NOTOK_MDF(fid=ncid)
      end if

      ! air mass
      if ( do_m ) then
        write (gol,'("  restore air mass ...")'); call goPr
        call MDF_Var_Par_Access( ncid, varid_m, MDF_INDEPENDENT, status )
        IF_NOTOK_MDF(fid=ncid)
        call MDF_Get_Var( ncid, varid_m, m_dat(n)%data(1:imr,1:jmr,1:lmr), status )
        IF_NOTOK_MDF(fid=ncid)
      end if

      !! surface fluxes
      !if ( do_sflux ) then
      !  write (gol,'("  restore surface fluxes ...")'); call goPr
      !  if ( slhf_dat(iglbsfc)%used ) then
      !    call MDF_Var_Par_Access( ncid, varid_slhf, MDF_INDEPENDENT, status )
      !    IF_NOTOK_MDF(fid=ncid)
      !    call MDF_Get_Var( ncid, varid_slhf, slhf_dat(iglbsfc)%data(1:nlon360,1:nlat180,1), status )
      !    IF_NOTOK_MDF(fid=ncid)
      !  end if
      !  if ( sshf_dat(iglbsfc)%used ) then
      !    call MDF_Var_Par_Access( ncid, varid_sshf, MDF_INDEPENDENT, status )
      !    IF_NOTOK_MDF(fid=ncid)
      !    call MDF_Get_Var( ncid, varid_sshf, sshf_dat(iglbsfc)%data(1:nlon360,1:nlat180,1), status )
      !    IF_NOTOK_MDF(fid=ncid)
      !  end if
      !end if

      ! tracer mass
      if ( do_rm ) then
      
        write (gol,'("  restore tracer mass ...")'); call goPr
        
        ! set independent data mode:
        call MDF_Var_Par_Access( ncid, varid_rm , MDF_INDEPENDENT, status )
        IF_NOTOK_MDF(fid=ncid)
#ifdef slopes   
        call MDF_Var_Par_Access( ncid, varid_rxm, MDF_INDEPENDENT, status )
        IF_NOTOK_MDF(fid=ncid)
        call MDF_Var_Par_Access( ncid, varid_rym, MDF_INDEPENDENT, status )
        IF_NOTOK_MDF(fid=ncid)
        call MDF_Var_Par_Access( ncid, varid_rzm, MDF_INDEPENDENT, status )
        IF_NOTOK_MDF(fid=ncid)
#endif
        if ( ntrace_chem > 0 ) then
           call MDF_Var_Par_Access( ncid, varid_rmc, MDF_INDEPENDENT, status )
           IF_NOTOK_MDF(fid=ncid)
        end if

        ! read list with tracer names in rcfile:        
        allocate( values_names(shp(2)) )
        call MDF_Get_Var( ncid, varid_names, values_names, status )
        IF_NOTOK_MDF(fid=ncid)

        ! loop over all tracers:
        do itr = 1, ntrace
       
          !WP! We skip this step as the tracer names in CTDAS are not unique
          !WP! Instead, simply assume the same number for each tracer in file and simulation

          ! search in list:
          !WP! call goMatchValue( names(itr), values_names, itr_file, status )

          itr_file=itr
          status = 0

          if ( status < 0 ) then  !WP! kind of obsolete now
             write(gol,'("Requested tracer `",a,"` not found in restart file!")') trim(names(itr)); call goErr
             IF_NOTOK_MDF(fid=ncid)
             !! TODO :
             !! insert (optional) default startup values here ??
             !! User can always add those in the restart file first,
             !! outside TM5... probably easier, faster, and more flexible.
          endif

          ! check ...
              if ( trim(names(itr)) /= trim(values_names(itr_file)) ) then
              write (gol,'("tracer `",a,"` does not have the same name as its counterpart `",a"` in the restart file")') trim(names(itr)), trim(values_names(itr_file)); call goErr
              status=1
              IF_NOTOK_MDF(fid=ncid)
          end if

          
          ! transported or short lived tracer ?
          if ( itr <= ntracet ) then

            ! check ...
             if ( itr_file > ntracet_restart ) then
              write (gol,'("tracer `",a,"` is transported but seems to be not-transported in restart file")') trim(names(itr)); call goErr
              status=1
              IF_NOTOK_MDF(fid=ncid)
            end if

            ! transported tracers distributed over tracers or levels:
            select case ( previous_par(n) )

              case ( 'tracer' )
        
                ! on this pe ?
                if ( tracer_active(itr) ) then
                
                  ! local tracer index:
                  itr_loc = tracer_loc(itr)

                  ! read tracer mass:
                  call MDF_Get_Var( ncid, varid_rm, &
                                mass_dat(n)%rm_t(1:imr,1:jmr,1:lmr,itr_loc:itr_loc), &
                                status, &
                                start=(/1,1,1,itr_file/), &
                                count=(/imr,jmr,lmr,1/) )
                  IF_NOTOK_MDF(fid=ncid)

#ifdef slopes
                  ! read slopes:
                  call MDF_Get_Var( ncid, varid_rxm, &
                                mass_dat(n)%rxm_t(1:imr,1:jmr,1:lmr,itr_loc:itr_loc), &
                                status, &
                                start=(/1,1,1,itr_file/), &
                                count=(/imr,jmr,lmr,1/) )
                  IF_NOTOK_MDF(fid=ncid)

                  call MDF_Get_Var( ncid, varid_rym, &
                                mass_dat(n)%rym_t(1:imr,1:jmr,1:lmr,itr_loc:itr_loc), &
                                status, &
                                start=(/1,1,1,itr_file/), &
                                count=(/imr,jmr,lmr,1/) )
                  IF_NOTOK_MDF(fid=ncid)

                  call MDF_Get_Var( ncid, varid_rzm, &
                                mass_dat(n)%rzm_t(1:imr,1:jmr,1:lmr,itr_loc:itr_loc), &
                                status, &
                                start=(/1,1,1,itr_file/), &
                                count=(/imr,jmr,lmr,1/) )
                  IF_NOTOK_MDF(fid=ncid)
#endif
                end if   ! tracer active

              case ( 'levels' )

                ! distributed over layers; might be empty on this pe
                if ( lmloc > 0 ) then

                  call MDF_Get_Var( ncid, varid_rm, &
                                  mass_dat(n)%rm_k(1:imr,1:jmr,1:lmloc,itr:itr), &
                                  status, &
                                  start=(/1,1,offsetl+1,itr_file/), &
                                  count=(/imr,jmr,lmloc,1/) )
                  IF_NOTOK_MDF(fid=ncid)

#ifdef slopes        
                  call MDF_Get_Var( ncid, varid_rxm, &
                                  mass_dat(n)%rxm_k(1:imr,1:jmr,1:lmloc,itr:itr), &
                                  status, &
                                  start=(/1,1,offsetl+1,itr_file/), &
                                  count=(/imr,jmr,lmloc,1/) )
                  IF_NOTOK_MDF(fid=ncid)

                  call MDF_Get_Var( ncid, varid_rym, &
                                  mass_dat(n)%rym_k(1:imr,1:jmr,1:lmloc,itr:itr), &
                                  status, &
                                  start=(/1,1,offsetl+1,itr_file/), &
                                  count=(/imr,jmr,lmloc,1/) )
                  IF_NOTOK_MDF(fid=ncid)

                  call MDF_Get_Var( ncid, varid_rzm, &
                                  mass_dat(n)%rzm_k(1:imr,1:jmr,1:lmloc,itr:itr), &
                                  status, &
                                  start=(/1,1,offsetl+1,itr_file/), &
                                  count=(/imr,jmr,lmloc,1/) )
                  IF_NOTOK_MDF(fid=ncid)
#endif
                end if   ! local layers > 0

              case default

                write (gol,'("unsupported distribution:")'); call goPr
                write (gol,'("  n         : ",i2)') n ;call goPr
                write (gol,'("  previous_par   : ",a)') previous_par(n) ;call goPr
                write (gol,'("in ",a)') rname; call goErr; status=1; return

            end select
          
          else   ! short lived tracer:
    
            ! check ...
             if ( itr_file <= ntracet_restart ) then
              write (gol,'("tracer `",a,"` is not-transported but seems to be transported in restart file")') trim(names(itr)); call goErr
              status=1
              IF_NOTOK_MDF(fid=ncid)
            end if

            ! non-transported tracer;
            ! distributed over layers, any layers on this pe ?
            if ( lmloc > 0 ) then

              ! read non-transported tracer into:
              !       chem_dat%rm_k(:,:,1:lmloc,ntracet+1:ntracet+ntrace_chem)
               call MDF_Get_Var( ncid, varid_rmc, &
                    chem_dat(n)%rm_k(1:imr,1:jmr,1:lmloc,itr:itr), &
                    status, &
                    start=(/1,1,offsetl+1,itr_file-ntracet_restart/), &
                              count=(/imr,jmr,lmloc,1/) )
              IF_NOTOK_MDF(fid=ncid)

            end if   ! lmloc > 0

          end if  ! transported or short-lived

        end do  ! tracers

        ! free mem for next region
        DeAllocate( values_names)      

     end if  ! do_rm

#ifdef with_tendencies
      ! production/loss/concentration
      if ( do_plc ) then
        write (gol,'(" restore tendency data ...")'); call goPr
        ! set independent data mode (not all processes are reading):
        do itr = 1, plc_ntr
          do ipr = 1, plc_npr
            call MDF_Var_Par_Access( ncid, varid_plc(itr,ipr), MDF_INDEPENDENT, status )
            IF_NOTOK_MDF(fid=ncid)
          end do
        end do
        ! distributed over levels:
        if ( lmloc > 0 ) then

          do itr = 1, plc_ntr
            do ipr = 1, plc_npr
          
              ! read plc slab:
              call MDF_Get_Var( ncid, varid_plc(itr,ipr), &
                              plc_dat(n,itr,ipr)%rm_k(1:imr,1:jmr,1:lmloc), status, &
                              start=(/1,1,offsetl+1/), count=(/imr,jmr,lmloc/) )
              IF_NOTOK_MDF(fid=ncid)

              ! get time attribute:
              call MDF_Get_Att( ncid, varid_plc(itr,ipr), 'time', time6, status )
              IF_NOTOK_MDF(fid=ncid)
              ! store time:
              call Set( plc_dat(n,itr,ipr)%t, time6=time6 )

            end do
          end do
          
        end if
        
      end if  ! do_plc
#endif

      ! o close file
      call MDF_Close( ncid, status )
      IF_NOTOK_RETURN(status=1)
    
    end do   ! regions
    
    ! ok
    status = 0
    
  end subroutine Restart_Read
  
  

end module Restart
