#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; if (isRoot) 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   ! read restart keys in rc file 
  public  ::  Restart_Done   ! nothing yet
  public  ::  Restart_Save   ! wrapper around Restart_Write
  public  ::  Restart_Write  ! write a restart file
  public  ::  Restart_Read   ! read a restart file
  public  ::  rs_write       ! model must write restart
  !
  ! !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. 
  !   10 Oct 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !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. Note that this does not solve all problems (but
  !   probably most of them): it is still possible that MDF_Close hangs...
  ! 
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RESTART_INIT
  !
  ! !DESCRIPTION: read settings from rcfile
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RESTART_INIT( status )
    !
    ! !USES:
    !
    use GO         , only : TrcFile, Init, Done, ReadRc
    use global_data, only : rcfile
    use global_data, only : outdir
    use meteodata  , only : lli
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)    ::  status
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = 'Restart_Init'
    type(TrcFile)                ::  rcF

    ! ---- begin
    
    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

    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

    status = 0

  END SUBROUTINE RESTART_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RESTART_DONE
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RESTART_DONE( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)    ::  status
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = 'Restart_Done'

    ! --- begin --------------------------------

    ! nothing to be done ...

    ! ok
    status = 0

  END SUBROUTINE RESTART_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  RESTART_SAVE
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RESTART_SAVE( status, extra, isfirst )
    !
    ! !USES:
    !
    use dims, only : idate
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)            ::  status
    !
    ! !INPUT PARAMETERS:
    !
    logical, intent(in), optional   ::  extra
    logical, intent(in), optional   ::  isfirst
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = 'Restart_Save'
    logical  :: is_extra
    real     :: t1, t2
    
    ! --- 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 cpu_time(t1)
          call Restart_Write( status, isfirst=isfirst )
          IF_NOTOK_RETURN(status=1)

          call cpu_time(t2)
          write (gol,*) " time to write restart [s]: ", t2-t1 ; call goPr

       end if  ! not extra

    end if  ! write at all

    ! ok
    status = 0

  END SUBROUTINE RESTART_SAVE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RESTART_FILENAME
  !
  ! !DESCRIPTION: Build restart filename from inputs.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RESTART_FILENAME( region, fname, status, key, dir, isfirst )
    !
    ! !USES:
    !
    use dims       , only : idate
    use global_data, only : outdir
    use meteodata  , only : lli
    !
    ! !INPUT PARAMETERS:
    !
    integer,          intent(in)            ::  region
    logical,          intent(in), optional  ::  isfirst
    character(len=*), intent(in), optional  ::  dir    
    character(len=*), intent(in), optional  ::  key
    !
    ! !OUTPUT PARAMETERS:
    !
    character(len=*), intent(out)   ::  fname
    integer,          intent(out)   ::  status
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
    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
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RESTART_WRITE
  !
  ! !DESCRIPTION: write restart
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RESTART_WRITE( status, key, region, isfirst )
    !
    ! !USES:
    !
    use GO          , only : Get
    use dims        , only : nregions, at, bt
    use chem_param  , only : ntracet, ntrace_chem, ntrace, names
    use partools    , only : isRoot
    use tm5_distgrid, only : dgrid,    Get_DistGrid, gather
    use global_data , only : mass_dat, chem_dat
    use meteodata   , only : global_lli, levi
    use meteodata   , only : sp_dat, phlb_dat, m_dat
    use MDF         , only : MDF_Create,  MDF_EndDef, MDF_Close
    use MDF         , only : MDF_Def_Dim, MDF_Def_Var
    use MDF         , only : MDF_Put_Att, MDF_Put_Var
    use MDF         , only : MDF_REPLACE, 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
    !   10 Oct 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = 'Restart_Write'

    integer               ::  imr, jmr, lmr, n
    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,       varid_at, varid_bt
    integer               ::  varid_sp,    varid_ph, varid_m
    integer               ::  varid_names, varid_rm
#ifdef slopes
    integer               ::  varid_rxm, varid_rym, varid_rzm
#endif
    integer               ::  varid_rmc
    integer               ::  rtype
    real, allocatable :: arr4d(:,:,:,:), arr3d(:,:,:)

    ! --- begin --------------------------------

    write (gol,'("write restart file(s) ...")'); call goPr

    ! loop over regions:
    REG: do n = 1, nregions

       ! only selected region ?
       if ( present(region) ) then
          if ( n /= region ) cycle
       end if

       ! entire region grid size
       imr = global_lli(n)%nlon
       jmr = global_lli(n)%nlat
       lmr = levi%nlev

       ! allocate 3D and 4D global arrays for gathering data
       if (isRoot) then
          allocate( arr4d(imr,jmr,lmr,ntracet) )
          allocate( arr3d(imr,jmr,lmr+1)       )
       else
          allocate( arr4d(1,1,1,1) )
          allocate( arr3d(1,1,1)   )
       endif

       ! name of restart file
       call Restart_FileName( n, fname, status, key=key, dir=rs_write_dir, isfirst=isfirst )
       IF_NOTOK_RETURN(status=1)

       write (gol,'("  destination : ",a)') trim(fname); call goPr

       if (isRoot) then
          !------------------
          ! OPEN NETCDF FILE
          !------------------        
          ! overwrite existing files (clobber)
          call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, ncid, status )
          IF_NOTOK_RETURN(status=1)

          !------------------
          ! DEFINE DIMENSIONS
          !------------------    
          call MDF_Def_Dim( ncid, 'lon', imr, dimid_lon, status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Def_Dim( ncid, 'lat', jmr, dimid_lat, status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Def_Dim( ncid, 'lev', lmr, dimid_lev, status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Def_Dim( ncid, 'hlev', lmr+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)

          !------------------
          ! DEFINE VARIABLES
          !------------------    
          select case ( kind(m_dat(n)%data) )
          case ( 4 ) ; rtype = MDF_FLOAT
          case ( 8 ) ; rtype = MDF_DOUBLE
          case default
             write (gol,'("unsupported real kind : ",i6)') kind(m_dat(n)%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

          ! 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

          !! 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

          !------------------
          ! END DEFINITION MODE
          !------------------    
          call MDF_EndDef( ncid, status )
          IF_NOTOK_MDF(fid=ncid)

       endif

       !------------------
       ! WRITE VARIABLES
       !------------------
       ! surface pressure
       call gather( dgrid(n), sp_dat(n)%data, arr3d(:,:,1:1), sp_dat(n)%halo, status)
       IF_NOTOK_RETURN(status=1)

       if (isRoot) call MDF_Put_Var( ncid, varid_sp, arr3d(:,:,1), status )
       IF_NOTOK_MDF(fid=ncid)

       ! half level pressure
       call gather( dgrid(n), phlb_dat(n)%data, arr3d, phlb_dat(n)%halo, status)
       IF_NOTOK_RETURN(status=1)

       if (isRoot) call MDF_Put_Var( ncid, varid_ph, arr3d, status)
       IF_NOTOK_MDF(fid=ncid)

       ! at, bt coefficients    
       if (isRoot) then
          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)
       end if

       ! air mass
       call gather( dgrid(n), m_dat(n)%data, arr4d(:,:,:,1), m_dat(n)%halo, status)
       IF_NOTOK_RETURN(status=1)

       if (isRoot) call MDF_Put_Var( ncid, varid_m, arr4d(:,:,:,1), 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:n360,1:n180,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:n360,1:n180,1), status )
       !IF_NOTOK_MDF(fid=ncid)

       ! tracer names
       if (isRoot) call MDF_Put_Var( ncid, varid_names, names, status )
       IF_NOTOK_MDF(fid=ncid)

       ! write transported tracers
       call gather( dgrid(n), mass_dat(n)%rm, arr4d, mass_dat(n)%halo, status)
       IF_NOTOK_RETURN(status=1)

       if (isRoot) call MDF_Put_Var( ncid, varid_rm, arr4d, status)
       IF_NOTOK_MDF(fid=ncid)

#ifdef slopes
       call gather( dgrid(n), mass_dat(n)%rxm, arr4d, mass_dat(n)%halo, status)
       IF_NOTOK_RETURN(status=1)

       if (isRoot) call MDF_Put_Var( ncid, varid_rxm, arr4d, status)
       IF_NOTOK_MDF(fid=ncid)

       call gather( dgrid(n), mass_dat(n)%rym, arr4d, mass_dat(n)%halo, status)
       IF_NOTOK_RETURN(status=1)

       if (isRoot) call MDF_Put_Var( ncid, varid_rym, arr4d, status)
       IF_NOTOK_MDF(fid=ncid)

       call gather( dgrid(n), mass_dat(n)%rzm, arr4d, mass_dat(n)%halo, status)
       IF_NOTOK_RETURN(status=1)

       if (isRoot) call MDF_Put_Var( ncid, varid_rzm, arr4d, status)
       IF_NOTOK_MDF(fid=ncid)
#endif

       ! write transported tracers
       call gather( dgrid(n), mass_dat(n)%rm, arr4d, mass_dat(n)%halo, status)
       IF_NOTOK_RETURN(status=1)

       if (isRoot) call MDF_Put_Var( ncid, varid_rm, arr4d, status)
       IF_NOTOK_MDF(fid=ncid)

#ifdef slopes
       call gather( dgrid(n), mass_dat(n)%rxm, arr4d, mass_dat(n)%halo, status)
       IF_NOTOK_RETURN(status=1)

       if (isRoot) call MDF_Put_Var( ncid, varid_rxm, arr4d, status)
       IF_NOTOK_MDF(fid=ncid)

       call gather( dgrid(n), mass_dat(n)%rym, arr4d, mass_dat(n)%halo, status)
       IF_NOTOK_RETURN(status=1)

       if (isRoot) call MDF_Put_Var( ncid, varid_rym, arr4d, status)
       IF_NOTOK_MDF(fid=ncid)

       call gather( dgrid(n), mass_dat(n)%rzm, arr4d, mass_dat(n)%halo, status)
       IF_NOTOK_RETURN(status=1)

       if (isRoot) call MDF_Put_Var( ncid, varid_rzm, arr4d, status)
       IF_NOTOK_MDF(fid=ncid)
#endif

       ! write non-transported tracers
       if (ntrace_chem > 0) then

          call gather( dgrid(n), chem_dat(n)%rm, arr4d(:,:,:,1:ntrace_chem),  chem_dat(n)%halo, status)
          IF_NOTOK_RETURN(status=1)

          if (isRoot) call MDF_Put_Var( ncid, varid_rmc, arr4d(:,:,:,1:ntrace_chem), status)
          IF_NOTOK_MDF(fid=ncid)

       end if

       ! Done
       if (isRoot) call MDF_Close( ncid, status )
       IF_NOTOK_RETURN(status=1)

       deallocate(arr4d, arr3d)

    end do REG

    status = 0

  END SUBROUTINE RESTART_WRITE
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RESTART_READ
  !
  ! !DESCRIPTION: Read restart file. Case of istart=33 (can read any of the
  !               available variables) or 32 (can read only tracer mass).
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RESTART_READ( status, region, &
                           surface_pressure, pressure,   air_mass,      surface_fluxes, &
                           tracer_mass,      tendencies, megan_history, nox_pulsing     )
    !
    ! !USES:
    !
    use GO,           only : TrcFile, Init, Done, ReadRc
    use GO,           only : goMatchValue
    use dims,         only : nregions, im, jm, istart, idate, idatei
    use grid,         only : TllGridInfo, TLevelInfo, Init, Done, Fill3D
    use chem_param,   only : ntracet, ntrace_chem, ntrace
    use chem_param,   only : names, tracer_name_len
    use partools,     only : isRoot, par_broadcast
    use tm5_distgrid, only : dgrid,  gather, scatter
    use global_data,  only : rcfile, mass_dat, chem_dat
    use meteodata,    only : levi, global_lli, 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,      MDF_COLLECTIVE
    use MDF,          only : MDF_Get_Att,        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
    !
    ! !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
    !   10 Oct 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:    
    !   - If we need to remap, then meteo is not read from restart.
    !     Airmass is still read but only to convert tracer masses to mixing ratios. 
    !     And istart should be 32.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Restart_Read'

    character(len=tracer_name_len), allocatable ::  values_names(:)
    character(len=256) ::  rs_read_dir, fname
    type(TrcFile)      ::  rcF    
    logical ::  exist
    logical ::  do_sp, do_ph, do_m, do_sflux, do_rm, do_megan, do_pulse
    integer ::  imr, jmr, lmr, imr_restart, jmr_restart, lmr_restart
    integer ::  n, region1, region2
    integer ::  ncid
    integer ::  varid_sp, varid_ph, varid_m, varid_rm, varid_rmc, varid_names
    !integer ::  varid_slhf, varid_sshf
    integer ::  itr, itr_file
    integer ::  ntracet_restart, dimid
    integer ::  shp(2)
#ifdef slopes
    integer ::  varid_rxm, varid_rym, varid_rzm
#endif
    ! global work arrays to read data
    real, allocatable :: tmp3d(:,:,:), airmass(:,:,:), run_airmass(:,:,:)
    real, allocatable :: rmt(:,:,:,:),rms(:,:,:,:), rmx(:,:,:,:),rmy(:,:,:,:), rmz(:,:,:,:)

    ! for remapping:
    logical :: need_vremap, need_hremap, need_remap
    integer :: varid_at, varid_bt
    real    :: dx, dy
    real, allocatable :: sp_gbl(:,:,:)
    real, allocatable :: at_restart(:), bt_restart(:)
    real, allocatable :: src_glb(:,:,:)
    type(TllGridInfo) :: lli_restart
    type(TLevelInfo)  :: levi_restart

    ! --- begin --------------------------------
    
    if ( istart /= 33 .and. istart /= 32 ) then
       write (gol,'("   skip read restart; istart not 33 or 32 but ",i2)') istart; call goPr
       status=0; return
    endif
    
    if ( any( idate /= idatei ) ) then
       write (gol,'("   skip read restart; idate not idatei but ",i4,5i2.2)') idate; call goPr
       status=0; return
    endif

    ! 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_rm      = .false.  ;  if ( present(tracer_mass      ) ) do_rm = tracer_mass
    do_m       = .false.  ;  if ( present(air_mass         ).and.(istart==33) ) do_m     = air_mass
    do_sp      = .false.  ;  if ( present(surface_pressure ).and.(istart==33) ) do_sp    = surface_pressure
    do_ph      = .false.  ;  if ( present(pressure         ).and.(istart==33) ) do_ph    = pressure
    do_sflux   = .false.  ;  if ( present(surface_fluxes   ).and.(istart==33) ) do_sflux = surface_fluxes
    do_megan   = .false.  ;  if ( present(megan_history    ).and.(istart==33) ) do_megan = megan_history
    do_pulse   = .false.  ;  if ( present(nox_pulsing      ).and.(istart==33) ) 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

    ! do we need anything?
    if(.not.(do_rm.or.do_m.or.do_sp.or.do_ph.or.do_sflux.or.do_megan.or.do_pulse))then
       status=0; return
    endif
    
    REG: do n = region1, region2

       imr = global_lli(n)%nlon
       jmr = global_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,'("   read restart file: ",a)') trim(fname); call goPr

       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

       ! ** open netcdf file
       if (isRoot) then 
          call MDF_Open( trim(fname), MDF_NETCDF4, MDF_READ, ncid, status )
          IF_NOTOK_RETURN(status=1)

          ! ** check for dimension compatibility
          call MDF_Inq_DimID( ncid, 'lev', dimid, status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Inquire_Dimension( ncid, dimid, status, length=lmr_restart )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Inq_DimID( ncid, 'lat', dimid, status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Inquire_Dimension( ncid, dimid, status, length=jmr_restart )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Inq_DimID( ncid, 'lon', dimid, status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Inquire_Dimension( ncid, dimid, status, length=imr_restart )
          IF_NOTOK_MDF(fid=ncid)
          need_vremap = (lmr /= lmr_restart)
          need_hremap = (jmr /= jmr_restart) .or. (imr /= imr_restart)
          need_remap  = need_hremap .or. need_vremap
       endif
       
       call par_broadcast( need_remap, status)
       IF_NOTOK_RETURN(status=1)
       
       if ((istart /= 32).and.need_remap) then
          status=1
          write(gol,*)' you must use istart=32 for using a restart file at different resolution'
          call goErr
          TRACEBACK; return
       endif

       ! work arrays
       if (isRoot) then
          allocate( rmt(imr,jmr,lmr,ntracet) )
          allocate( rmx(imr,jmr,lmr,ntracet) )
          allocate( rmy(imr,jmr,lmr,ntracet) )
          allocate( rmz(imr,jmr,lmr,ntracet) )
          if ( ntrace_chem > 0 ) allocate( rms(imr,jmr,lmr,ntracet+1:ntracet+ntrace_chem) )
          allocate( tmp3d(imr,jmr,lmr+1    ) )
          allocate( airmass(imr_restart,jmr_restart,lmr_restart    ) )
          if (istart==32) allocate( run_airmass(imr,jmr,lmr) )
       else
          allocate( rmt(1,1,1,1) )
          allocate( rmx(1,1,1,1) )
          allocate( rmy(1,1,1,1) )
          allocate( rmz(1,1,1,1) )
          if ( ntrace_chem > 0 ) allocate( rms(1,1,1,1) )
          allocate( airmass(1,1,1) )
          if (istart==32) allocate( run_airmass(1,1,1) )
          allocate( tmp3d(1,1,1) )
       endif
       
       if (istart==32) then
         CALL GATHER( dgrid(n), m_dat(n)%data, run_airmass, m_dat(n)%halo, status )
         IF_NOTOK_RETURN(status=1)
       endif
       
       ! prepare for remap 
       if (need_remap .and. do_rm) then

          write (gol,'("    remap tracer from restart file")') ; call goPr
          
          if (isRoot) then
             allocate( sp_gbl(imr,jmr,1) )
             allocate( src_glb(imr_restart,jmr_restart,lmr_restart))
          else
             allocate(sp_gbl(1,1,1))
             allocate(src_glb(1,1,1))
          endif

          call gather( dgrid(n), sp_dat(n)%data, sp_gbl, sp_dat(n)%halo, status)
          IF_NOTOK_RETURN(status=1)

          ! init to 0 in case of data not found in file
          rmt=0.
          rms=0.

          ! init lli_restart/levi_restart
          if (isRoot) then
            allocate(at_restart(lmr_restart+1))
            allocate(bt_restart(lmr_restart+1))
            !
            call MDF_Inq_VarID( ncid, 'at', varid_at, status )
            IF_NOTOK_MDF(fid=ncid)
            !
            call MDF_Get_Var( ncid, varid_at, at_restart(1:(lmr_restart+1)), status )
            IF_NOTOK_MDF(fid=ncid)
            !
            call MDF_Inq_VarID( ncid, 'bt', varid_bt, status )
            IF_NOTOK_MDF(fid=ncid)
            !
            call MDF_Get_Var( ncid, varid_bt, bt_restart(1:(lmr_restart+1)), status )
            IF_NOTOK_MDF(fid=ncid)
            !
            call Init( levi_restart, lmr_restart, at_restart, bt_restart, status )
            IF_NOTOK_RETURN(status=1)
            !
            deallocate(at_restart,bt_restart)
            !
            dx=360./imr_restart
            dy=180./jmr_restart
            call Init( lli_restart, -180.+0.5*dx, dx, imr_restart, &
                                     -90.+0.5*dy, dy, jmr_restart, status )
            IF_NOTOK_RETURN(status=1)
          endif
       endif
              
       ! ** get variables id
       if (isRoot) then
          ! surface pressure
          if ( do_sp ) call MDF_Inq_VarID( ncid, 'sp', varid_sp, status )
          IF_NOTOK_MDF(fid=ncid)

          ! half level pressure
          if ( do_ph ) call MDF_Inq_VarID( ncid, 'ph', varid_ph, status )
          IF_NOTOK_MDF(fid=ncid)

          ! air mass
          call MDF_Inq_VarID( ncid, 'm', varid_m, status )
          IF_NOTOK_MDF(fid=ncid)

          !! surface fluxes
          !if ( do_sflux ) then
          !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)

             ! tracers mass 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 if any
             if ( ntrace_chem > 0 ) then
                call MDF_Inq_VarID( ncid, 'rmc', varid_rmc, status )
                IF_NOTOK_MDF(fid=ncid)                    
             end if
          end if
       end if

       ! *** READ VARIABLES ***
       
       if ( do_sp ) then
          write (gol,'("    restore surface pressure ...")'); call goPr

          if (isRoot) call MDF_Get_Var( ncid, varid_sp, tmp3d(:,:,1), status )
          IF_NOTOK_MDF(fid=ncid)

          call scatter( dgrid(n), sp_dat(n)%data, tmp3d(:,:,1:1), sp_dat(n)%halo, status)
          IF_NOTOK_RETURN(status=1)

       end if

       if ( do_ph ) then
          write (gol,'("    restore half level pressure ...")'); call goPr

          if (isRoot) call MDF_Get_Var( ncid, varid_ph, tmp3d, status)
          IF_NOTOK_MDF(fid=ncid)

          call scatter( dgrid(n), phlb_dat(n)%data, tmp3d, phlb_dat(n)%halo, status)
          IF_NOTOK_RETURN(status=1)
       end if
       
       ! get air mass in all cases
       if (isRoot) call MDF_Get_Var( ncid, varid_m, airmass, status)
       IF_NOTOK_MDF(fid=ncid)

       if ( do_m ) then
          write (gol,'("    restore air mass ...")'); call goPr
         
          call scatter( dgrid(n), m_dat(n)%data, airmass, m_dat(n)%halo, status)
          IF_NOTOK_RETURN(status=1)
       end if


       ! tracer mass
       READRM: if ( do_rm ) then

          write (gol,'("    restore tracer mass ...")'); call goPr

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

          ! loop over all model tracers
          do itr = 1, ntrace

             if (isRoot) then

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

                if ( status < 0 ) then
                   write(gol,'("*WARNING* Requested tracer `",a,"` not FOUND in restart file!")') trim(names(itr))

                   if (istart /= 32) then
                      call goErr
                      IF_NOTOK_MDF(fid=ncid)
                   else
                      status=0
                      call goPr
                      if ( itr <= ntracet ) then
                         rmt(:,:,:,itr) = 1.e-30
                         write(gol,'("*WARNING* Requested TRANSPORTED tracer `",a,"` has been SET to a default value of 1.e-30")') trim(names(itr))
                      else
                         rms(:,:,:,itr) = 1.e-30
                         write(gol,'("*WARNING* Requested SHORT-LIVED tracer `",a,"` has been SET to a default value of 1.e-30")') trim(names(itr))
                      endif
                      call goPr
                   endif
                else

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

                      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

                      if (need_remap) then
                         call MDF_Get_Var( ncid, varid_rm, src_glb, status, start=(/1,1,1,itr_file/))
                         IF_NOTOK_MDF(fid=ncid)
                         
                         src_glb = src_glb / airmass
                         
                         call Fill3D( global_lli(n), levi, 'n', sp_gbl(:,:,1), rmt(:,:,:,itr), &
                                      lli_restart,   levi_restart, src_glb, 'mass-aver', status )
                         IF_NOTOK_RETURN(status=1)
                         
                         rmt(:,:,:,itr) = rmt(:,:,:,itr) * run_airmass
                         
                      else
                         call MDF_Get_Var( ncid, varid_rm, rmt(:,:,:,itr), status, start=(/1,1,1,itr_file/))
                         IF_NOTOK_MDF(fid=ncid)

                         if (istart==32) then
                           rmt(:,:,:,itr) = rmt(:,:,:,itr) * run_airmass / airmass
                         endif
                         
                      endif
#ifdef slopes
                      ! read slopes
                      if ((.not. need_remap) .and. (istart==33)) then

                         if (isRoot) call MDF_Get_Var( ncid, varid_rxm, rmx(:,:,:,itr), status, start=(/1,1,1,itr_file/))
                         IF_NOTOK_MDF(fid=ncid)                   

                         if (isRoot) call MDF_Get_Var( ncid, varid_rym, rmy(:,:,:,itr), status, start=(/1,1,1,itr_file/))
                         IF_NOTOK_MDF(fid=ncid)                   

                         if (isRoot) call MDF_Get_Var( ncid, varid_rzm, rmz(:,:,:,itr), status, start=(/1,1,1,itr_file/))
                         IF_NOTOK_MDF(fid=ncid)

                      endif
#endif
                   else   ! short lived tracer:

                      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

                      itr_file = itr_file - ntracet_restart
                      if (need_remap) then
                         call MDF_Get_Var( ncid, varid_rmc, src_glb, status, start=(/1,1,1,itr_file/) )
                         IF_NOTOK_MDF(fid=ncid)
                         
                         src_glb = src_glb / airmass
                         
                         call Fill3D( global_lli(n), levi, 'n', sp_gbl(:,:,1), rms(:,:,:,itr), &
                                      lli_restart,   levi_restart, src_glb, 'mass-aver', status )
                         IF_NOTOK_RETURN(status=1)

                         rms(:,:,:,itr) = rms(:,:,:,itr) * run_airmass
                         
                      else
                         call MDF_Get_Var( ncid, varid_rmc, rms(:,:,:,itr), status, start=(/1,1,1,itr_file/) )
                         IF_NOTOK_MDF(fid=ncid)

                         if (istart==32) then
                           rms(:,:,:,itr) = rms(:,:,:,itr) * run_airmass / airmass
                         endif

                      endif

                   end if  ! transported or short-lived
                endif ! in the file
             endif ! root
          end do  ! tracers

          ! distribute
          call scatter( dgrid(n), mass_dat(n)%rm, rmt, mass_dat(n)%halo, status)
          IF_NOTOK_RETURN(status=1)

          if ( ntrace_chem > 0 ) then
             call scatter( dgrid(n), chem_dat(n)%rm, rms, chem_dat(n)%halo, status)
             IF_NOTOK_RETURN(status=1)
          endif

#ifdef slopes
          if ((.not. need_remap).and.(istart==33)) then
             call scatter( dgrid(n), mass_dat(n)%rxm, rmx, mass_dat(n)%halo, status)
             IF_NOTOK_RETURN(status=1)
             call scatter( dgrid(n), mass_dat(n)%rym, rmy, mass_dat(n)%halo, status)
             IF_NOTOK_RETURN(status=1)
             call scatter( dgrid(n), mass_dat(n)%rzm, rmz, mass_dat(n)%halo, status)
             IF_NOTOK_RETURN(status=1)
          else
             ! Ensure that slopes are initialized to "unset" values of 0.0.  Wouter says that
             ! we could remap levels for rxm et al., but 0.0 will also work.  The noise
             ! induced from remapping the rm array is almost certainly bigger than any issues
             ! from having this "default=0.0" slopes information.  -ARJ 1 Jan 12
             mass_dat(n)%rxm = 0.0                     
             mass_dat(n)%rym = 0.0                     
             mass_dat(n)%rzm = 0.0                     
          endif
#endif
          ! free mem for next region
          deallocate( values_names)      

          if (need_remap) then
            deallocate(sp_gbl,src_glb)
            if (isRoot) then
              call Done( levi_restart, status )
              IF_NOTOK_RETURN(status=1)
              call Done( lli_restart, status )
              IF_NOTOK_RETURN(status=1)
            endif
          endif
          
       ENDIF READRM

       ! clean
       deallocate(rmt)
       if ( ntrace_chem > 0 ) deallocate(rms)
#ifdef slopes
       deallocate(rmx, rmy, rmz)
#endif       
       deallocate( tmp3d )
       deallocate( airmass)
       if (istart==32) deallocate(run_airmass)

       if (isRoot) call MDF_Close( ncid, status )
       IF_NOTOK_RETURN(status=1)

    ENDDO REG

    status = 0

  END SUBROUTINE RESTART_READ
  !EOC

END MODULE RESTART