ecearth3/sources/tm5mp/proj/cbm4/tm5_restart.F90

!
#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. This version differs from the
!               'base' version by accounting for "with_online_nox", 
!               "with_online_bvoc", and "with_m7" cpp flags, which read/write
!               additional datasets.
!\\
!\\
! !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: 
  !   25 Aug 2010 - P. Le Sager - Merged with Base version for Pycasso 
  !    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: 
    !   24 Aug 2010 - P. Le Sager - merged w/ trunk for pycasso
    !
    !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 dims        , only : iglbsfc
    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
#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
#ifdef with_online_bvoc
    use emission_bvoc_data, only : megan, pceea
    use emission_bvoc_data, only : ndays_history, nhours_history, n_layers
    use emission_bvoc_data, only : skt_daily, pdir_daily, pdif_daily, ssr_daily
    use emission_bvoc_data, only : skt_10d_history, pdir_10d_history, pdif_10d_history, ssr_10d_history
    use emission_bvoc_data, only : skt_hourly, pdir_hourly, pdif_hourly
    use emission_bvoc_data, only : skt_24h_history, pdir_24h_history, pdif_24h_history
    use chem_param,         only : iisop
#endif
#ifdef with_online_nox
    use online_nox_data, only : pulsing_on, ndrydays
    use online_nox_data, only : cp_daily, lsp_daily
    use online_nox_data, only : cp_history, lsp_history
    use online_nox_data, only : pulsing_field, plsday_field, plsdurat_field
#endif
#ifdef with_m7 
    use mo_aero_m7,          only : nsol, nmod
    use m7_data,             only : h2o_mode, rw_mode, rwd_mode
#endif
    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
    !   16 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !   - Serial writing not tested for cases: with_online_nox and with_online_bvoc (February 19, 2018)
    !
    !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
#ifdef with_tendencies
    integer               ::  varid_plc(plc_ntr,plc_npr)
    integer               ::  itr, ipr
    integer               ::  time6(6)
#endif
#ifdef with_online_bvoc 
    integer               ::  dimid_lon_bvoc, dimid_lat_bvoc
    integer               ::  dimid_days_history, dimid_hours_history, dimid_layers
    integer               ::  varid_skt_daily, varid_pdir_daily, varid_pdif_daily, varid_ssr_daily
    integer               ::  varid_skt_10d_history, varid_pdir_10d_history, varid_pdif_10d_history, varid_ssr_10d_history
    integer               ::  varid_skt_hourly, varid_pdir_hourly, varid_pdif_hourly
    integer               ::  varid_skt_24h_history, varid_pdir_24h_history, varid_pdif_24h_history
#endif
#ifdef with_online_nox
    integer               ::  dimid_lon_nox, dimid_lat_nox
    integer               ::  dimid_drydays
    integer               ::  varid_cp_daily, varid_lsp_daily
    integer               ::  varid_cp_history, varid_lsp_history
    integer               ::  varid_pulsing, varid_plsday, varid_plsdurat
#endif
#ifdef with_m7
    integer               ::  varid_h2o, varid_rw, varid_rwd
    integer               ::  dimid_nsol, dimid_nmod
    integer               ::  imode
    character(len=3), parameter :: h2o_name = 'h2o'
    character(len=3), parameter :: rwd_name = 'rwd'
    character(len=2), parameter :: rw_name  = 'rw'
#endif
    integer               ::  rtype, n360, n180

    real, allocatable :: arr4d(:,:,:,:), arr3d(:,:,:)
#if defined(with_online_bvoc) || defined(with_online_nox)
    real, allocatable :: glb_sfc3D(:,:,:), glb_sfc4D(:,:,:,:)
#endif  

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

       ! get extra bounds for 1x1 dataset
#if defined(with_online_bvoc) || defined(with_online_nox)
       if(n==1) then
          n360 = dgrid(iglbsfc)%im_region
          n180 = dgrid(iglbsfc)%jm_region
          if (isRoot) then
             allocate( glb_sfc3d(n360, n180, max(ndays_history, n_layers, ndrydays, nhours_history ) )
             allocate( glb_sfc4d(n360, n180, n_layers, ndays_history)    )
          else
             allocate( glb_sfc3d(1,1,1)   )
             allocate( glb_sfc4d(1,1,1,1) )
          endif
       end if
#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, 'lon_sfc', global_lli(iglbsfc)%nlon, dimid_lon_sfc, status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Def_Dim( ncid, 'lat_sfc', global_lli(iglbsfc)%nlat, dimid_lat_sfc, 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)

#ifdef with_online_bvoc
          ! MEGAN/PCEEA history
          if ( (megan .or. pceea) .and. (n == 1) ) then
             call MDF_Def_Dim( ncid, 'lon_bvoc', n360, dimid_lon_bvoc ,status)
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Def_Dim( ncid, 'lat_bvoc', n180, dimid_lat_bvoc ,status)
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Def_Dim( ncid, 'day_bvoc', ndays_history, dimid_days_history ,status)
             IF_NOTOK_MDF(fid=ncid)
          endif
          if (megan .and. (n == 1) ) then
             call MDF_Def_Dim( ncid, 'hour_bvoc', nhours_history, dimid_hours_history ,status)
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Def_Dim( ncid, 'layer_bvoc', n_layers, dimid_layers ,status)
             IF_NOTOK_MDF(fid=ncid)
          endif
#endif 

#ifdef with_online_nox
          ! precipitation history and pulsing parameters
          if (pulsing_on .and. (n == 1) ) then
             call MDF_Def_Dim( ncid, 'lon_nox', n360, dimid_lon_nox ,status)
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Def_Dim( ncid, 'lat_nox', n180, dimid_lat_nox ,status)
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Def_Dim( ncid, 'day_nox', ndrydays, dimid_drydays ,status)
             IF_NOTOK_MDF(fid=ncid)
          endif
#endif

#ifdef with_m7
          ! --------------------
          ! M7 fields for optics
          ! --------------------
          call MDF_Def_Dim( ncid, 'nsol', nsol, dimid_nsol, status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Def_Dim( ncid, 'nmod', nmod, dimid_nmod, status )
          IF_NOTOK_MDF(fid=ncid)
#endif

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

#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

#ifdef with_online_bvoc
          ! MEGAN/PCEEA history parameters
          if ( (megan .or. pceea) .and. (n == 1) ) then
             call MDF_Def_Var( ncid, 'skt_daily', rtype, (/dimid_lon_bvoc,dimid_lat_bvoc/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', 'average skin temperature since the start of this day', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'K', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_skt_daily = varid
             !
             call MDF_Def_Var( ncid, 'skt_10d_history', rtype, (/dimid_lon_bvoc,dimid_lat_bvoc,dimid_days_history/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', '10-day skin temperature record', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'K', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_skt_10d_history = varid
          endif
          if (megan .and. (n == 1) ) then
             !
             call MDF_Def_Var( ncid, 'pdir_daily', rtype, (/dimid_lon_bvoc,dimid_lat_bvoc/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', &
                  'average direct component of the photosynthetic photon flux density since the start of this day', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'umol/(m2 s)', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_pdir_daily = varid
             !
             call MDF_Def_Var( ncid, 'pdir_10d_history', rtype, (/dimid_lon_bvoc,dimid_lat_bvoc,dimid_days_history/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', &
                  '10-day record of the direct component of the photosynthetic photon flux density', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'umol/(m2 s)', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_pdir_10d_history = varid
             !
             call MDF_Def_Var( ncid, 'pdif_daily', rtype, (/dimid_lon_bvoc,dimid_lat_bvoc,dimid_layers/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', &
                  'average diffuse component of the photosynthetic photon flux density since the start of this day', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'umol/(m2 s)', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_pdif_daily = varid
             !
             call MDF_Def_Var( ncid, 'pdif_10d_history', rtype, &
                  (/dimid_lon_bvoc,dimid_lat_bvoc,dimid_layers,dimid_days_history/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', &
                  '10-day record of the diffuse component of the photosynthetic photon flux density', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'umol/(m2 s)', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_pdif_10d_history = varid
             !
             call MDF_Def_Var( ncid, 'skt_hourly', rtype, (/dimid_lon_bvoc,dimid_lat_bvoc/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', 'averaged skin temperature since the start of this hour', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'K', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_skt_hourly = varid
             !
             call MDF_Def_Var( ncid, 'skt_24h_history', rtype, (/dimid_lon_bvoc,dimid_lat_bvoc,dimid_hours_history/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', '24-hour skin temperature record', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'K', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_skt_24h_history = varid
             !
             call MDF_Def_Var( ncid, 'pdir_hourly', rtype, (/dimid_lon_bvoc,dimid_lat_bvoc/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', &
                  'average direct component of the photosynthetic photon flux density since the start of this hour', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'umol/(m2 s)', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_pdir_hourly = varid
             !
             call MDF_Def_Var( ncid, 'pdir_24h_history', rtype, (/dimid_lon_bvoc,dimid_lat_bvoc,dimid_hours_history/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', &
                  '24-hour record of the direct component of the photosynthetic photon flux density', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'umol/(m2 s)', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_pdir_24h_history = varid
             !
             call MDF_Def_Var( ncid, 'pdif_hourly', rtype, (/dimid_lon_bvoc,dimid_lat_bvoc,dimid_layers/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', &
                  'average diffuse component of the photosynthetic photon flux density since the start of this hour', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'umol/(m2 s)', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_pdif_hourly = varid
             !
             call MDF_Def_Var( ncid, 'pdif_24h_history', rtype, &
                  (/dimid_lon_bvoc,dimid_lat_bvoc,dimid_layers,dimid_hours_history/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', &
                  '24-hour record of the diffuse component of the photosynthetic photon flux density', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'umol/(m2 s)', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_pdif_24h_history = varid
             !
          else if ( pceea .and. (n == 1) ) then 
             !
             call MDF_Def_Var( ncid, 'ssr_daily', rtype, (/dimid_lon_bvoc,dimid_lat_bvoc/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', 'average surface solar radiation since the start of this day', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'W/m2', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_ssr_daily = varid
             !
             call MDF_Def_Var( ncid, 'ssr_10d_history', rtype, (/dimid_lon_bvoc,dimid_lat_bvoc,dimid_days_history/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', '10-day surface solar radiation record', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'W/m2', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_ssr_10d_history = varid
             !
          endif

#endif 

#ifdef with_online_nox
          ! precipitation history and pulsing parameters
          if (pulsing_on .and. (n == 1) ) then
             !
             call MDF_Def_Var( ncid, 'cp_daily', rtype, (/dimid_lon_nox,dimid_lat_nox/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', 'daily accumulated convective rainfall', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'm', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_cp_daily = varid
             !
             call MDF_Def_Var( ncid, 'lsp_daily', rtype, (/dimid_lon_nox,dimid_lat_nox/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', 'daily accumulated large-scale rainfall', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'm', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_lsp_daily = varid
             !
             call MDF_Def_Var( ncid, 'cp_history', rtype, (/dimid_lon_nox,dimid_lat_nox,dimid_drydays/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', '14-day convective rainfall record', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'm', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_cp_history = varid
             !
             call MDF_Def_Var( ncid, 'lsp_history', rtype, (/dimid_lon_nox,dimid_lat_nox,dimid_drydays/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', '14-day large-scale rainfall record', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'm', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_lsp_history = varid
             !
             call MDF_Def_Var( ncid, 'pulsing', rtype, (/dimid_lon_nox,dimid_lat_nox/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', 'pulsing regime', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'unity', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_pulsing = varid
             !
             call MDF_Def_Var( ncid, 'plsday', rtype, (/dimid_lon_nox,dimid_lat_nox/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', 'time of pulsing', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'days', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_plsday = varid
             !
             call MDF_Def_Var( ncid, 'plsdurat', rtype, (/dimid_lon_nox,dimid_lat_nox/), varid, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'long_name', 'duration of pulse', status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Put_Att( ncid, varid, 'unit', 'days', status )
             IF_NOTOK_MDF(fid=ncid)
             varid_plsdurat = varid
             !
          endif
#endif
#ifdef with_m7
#ifndef without_chemistry
          ! --------------------
          ! M7 fields for optics
          ! --------------------

          ! water fields
          call MDF_Def_Var( ncid, trim(h2o_name), rtype, (/dimid_lon,dimid_lat,dimid_lev,dimid_nsol/), varid, status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Put_Att( ncid, varid, 'long_name', 'aerosol water content', status)
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Put_Att( ncid, varid, 'unit', 'kg', status)
          IF_NOTOK_MDF(fid=ncid)
          varid_h2o = varid

          ! dry radii for soluble modes
          call MDF_Def_Var( ncid, trim(rwd_name), rtype, (/dimid_lon,dimid_lat,dimid_lev,dimid_nsol/), varid, status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Put_Att( ncid, varid, 'long_name', 'mode dry radius', status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Put_Att( ncid, varid, 'unit', 'm', status )
          IF_NOTOK_MDF(fid=ncid)
          varid_rwd = varid

          ! wet radii
          call MDF_Def_Var( ncid, trim(rw_name), rtype, (/dimid_lon,dimid_lat,dimid_lev,dimid_nmod/), varid, status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Put_Att( ncid, varid, 'long_name', 'mode radius', status )
          IF_NOTOK_MDF(fid=ncid)
          call MDF_Put_Att( ncid, varid, 'unit', 'm', status )
          IF_NOTOK_MDF(fid=ncid)
          varid_rw = varid
#endif
#endif
          !------------------
          ! 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 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

#ifdef with_tendencies

       ! write production/loss/concentration levels on this pe:
       do itr = 1, plc_ntr
          do ipr = 1, plc_npr

             call gather( dgrid(n), plc_dat(n,itr,ipr)%rm, arr3d(:,:,1:lmr), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_plc(itr,ipr), arr3d(:,:,1:lmr), status )
             IF_NOTOK_MDF(fid=ncid)

          end do
       end do
#endif

#ifdef with_online_bvoc
       ! MEGAN/PCEEA history parameters; write only once
       if (n==1) then
          if (megan .or. pceea) then

             call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_skt_daily, glb_sfc3D(:,:,1), status)
             IF_NOTOK_MDF(fid=ncid)

             call gather( dgrid(n), skt_10d_history, glb_sfc3D(:,:,1:ndays_history), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_skt_10d_history, glb_sfc3D(:,:,1:ndays_history), status)
             IF_NOTOK_MDF(fid=ncid)

          endif
          !
          if (megan) then

             call gather( dgrid(n), pdir_daily, glb_sfc3D(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_pdir_daily, glb_sfc3D(:,:,1), status)
             IF_NOTOK_MDF(fid=ncid)

             call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1:ndays_history), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_pdir_10d_history, glb_sfc3D(:,:,1:ndays_history), status)
             IF_NOTOK_MDF(fid=ncid)

             call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1:n_layers), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_pdif_daily, glb_sfc3D(:,:,1:n_layers), status)
             IF_NOTOK_MDF(fid=ncid)

             call gather( dgrid(n), skt_daily, glb_sfc4D, 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_pdif_10d_history, glb_sfc4D, status)
             IF_NOTOK_MDF(fid=ncid)

             call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_skt_hourly, glb_sfc3D(:,:,1), status)
             IF_NOTOK_MDF(fid=ncid)

             call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1:nhours_history), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_skt_24h_history,  glb_sfc3D(:,:,1:nhours_history), status)
             IF_NOTOK_MDF(fid=ncid)

             call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_pdir_hourly, glb_sfc3D(:,:,1), status)
             IF_NOTOK_MDF(fid=ncid)

             call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1:nhours_history), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_pdir_24h_history, glb_sfc3D(:,:,1:nhours_history), status)
             IF_NOTOK_MDF(fid=ncid)

             call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1:n_layers), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_pdif_hourly, glb_sfc3D(:,:,1:n_layers), status)
             IF_NOTOK_MDF(fid=ncid)

             call gather( dgrid(n), skt_daily, glb_sfc4D, 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_pdif_24h_history, glb_sfc4D, status)
             IF_NOTOK_MDF(fid=ncid)

          else if (pceea) then

             call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_ssr_daily, glb_sfc3D(:,:,1), status)
             IF_NOTOK_MDF(fid=ncid)

             call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1:ndays_history), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Put_Var( ncid, varid_ssr_10d_history, glb_sfc3D(:,:,1:ndays_history), status)
             IF_NOTOK_MDF(fid=ncid)

          endif
          ! 
       endif
#endif

#ifdef with_online_nox
       if (pulsing_on .and. (n == 1) ) then
          ! precipitation history and pulsing parameters; write only once

          call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1), 0, status)
          IF_NOTOK_RETURN(status=1)

          if (isRoot) call MDF_Put_Var( ncid,  varid_cp_daily, glb_sfc3D(:,:,1), status)
          IF_NOTOK_MDF(fid=ncid)

          call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1), 0, status)
          IF_NOTOK_RETURN(status=1)

          if (isRoot) call MDF_Put_Var( ncid,  varid_lsp_daily,   glb_sfc3D(:,:,1), status)
          IF_NOTOK_MDF(fid=ncid)

          call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1:ndrydays), 0, status)
          IF_NOTOK_RETURN(status=1)

          if (isRoot) call MDF_Put_Var( ncid,  varid_cp_history,  glb_sfc3D(:,:,1:ndrydays),  status)
          IF_NOTOK_MDF(fid=ncid)

          call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1:ndrydays), 0, status)
          IF_NOTOK_RETURN(status=1)

          if (isRoot) call MDF_Put_Var( ncid,  varid_lsp_history, glb_sfc3D(:,:,1:ndrydays), status)
          IF_NOTOK_MDF(fid=ncid)

          call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1), 0, status)
          IF_NOTOK_RETURN(status=1)

          if (isRoot) call MDF_Put_Var( ncid,  varid_pulsing, glb_sfc3D(:,:,1), status)
          IF_NOTOK_MDF(fid=ncid)

          call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1), 0, status)
          IF_NOTOK_RETURN(status=1)

          if (isRoot) call MDF_Put_Var( ncid,  varid_plsday, glb_sfc3D(:,:,1), status)
          IF_NOTOK_MDF(fid=ncid)

          call gather( dgrid(n), skt_daily, glb_sfc3D(:,:,1), 0, status)
          IF_NOTOK_RETURN(status=1)

          if (isRoot) call MDF_Put_Var( ncid,  varid_plsdurat, glb_sfc3D(:,:,1), status)
          IF_NOTOK_MDF(fid=ncid)

       endif
#endif

#ifdef with_m7
#ifndef without_chemistry        

       do imode=1,nsol   
          call gather( dgrid(n), h2o_mode(n,imode)%d3, arr4d(:,:,:,imode), h2o_mode(n,imode)%halo, status)
          IF_NOTOK_RETURN(status=1)
       enddo

       if (isRoot) call MDF_Put_Var( ncid, varid_h2o, arr4d(:,:,:,1:nsol), status)
       IF_NOTOK_MDF(fid=ncid)

       do imode=1,nsol   
          call gather( dgrid(n), rwd_mode(n,imode)%d3, arr4d(:,:,:,imode), rwd_mode(n,imode)%halo, status)
          IF_NOTOK_RETURN(status=1)
       enddo

       if (isRoot) call MDF_Put_Var( ncid, varid_rwd, arr4d(:,:,:,1:nsol), status)
       IF_NOTOK_MDF(fid=ncid)

       do imode=1,nmod
          call gather( dgrid(n), rw_mode(n,imode)%d3, arr4d(:,:,:,imode), rw_mode(n,imode)%halo, status)
          IF_NOTOK_RETURN(status=1)
       enddo

       if (isRoot) call MDF_Put_Var( ncid, varid_rw, arr4d(:,:,:,1:nmod), status)
       IF_NOTOK_MDF(fid=ncid)
#endif
#endif

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

       deallocate(arr4d, arr3d)

#if defined(with_online_bvoc) || defined(with_online_nox)
       deallocate(glb_sfc3D, glb_sfc4D)
#endif  

    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 dims,         only : iglbsfc
    use grid,         only : TllGridInfo, TLevelInfo, Init, Done, Fill3D
    use chem_param,   only : ntracet, ntrace_chem, ntrace, ich4
    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
#ifdef with_online_bvoc
    use emission_bvoc_data, only : megan, pceea
    use emission_bvoc_data, only : ndays_history, nhours_history, n_layers
    use emission_bvoc_data, only : skt_daily, pdir_daily, pdif_daily, ssr_daily
    use emission_bvoc_data, only : skt_10d_history, pdir_10d_history, pdif_10d_history, ssr_10d_history
    use emission_bvoc_data, only : skt_hourly, pdir_hourly, pdif_hourly
    use emission_bvoc_data, only : skt_24h_history, pdir_24h_history, pdif_24h_history
    use chem_param,         only : iisop
    use partools,           only : tracer_active
#endif
#ifdef with_online_nox
    use online_nox_data, only : ndrydays
    use online_nox_data, only : cp_daily, lsp_daily
    use online_nox_data, only : cp_history, lsp_history
    use online_nox_data, only : pulsing_field, plsday_field, plsdurat_field
!    use partools,        only : root_k
#endif
#ifdef with_m7
    use mo_aero_m7,      only : nsol, nmod 
    use m7_data,         only : h2o_mode, rw_mode, rwd_mode
#endif 
    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
    !   16 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !   - logically, if we need to remap, then meteo is not read from restart
    !      (but from met field and used for remapping): in other words, only
    !      tracer mass is read, and istart should be 32.
    !   - Serial reading not tested for cases: with_online_nox and with_online_bvoc (February 19, 2018)
    !
    !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
    
#ifdef with_online_bvoc
    integer               ::  varid_skt_daily, varid_pdir_daily, varid_pdif_daily, varid_ssr_daily
    integer               ::  varid_skt_10d_history, varid_pdir_10d_history, varid_pdif_10d_history, varid_ssr_10d_history
    integer               ::  varid_skt_hourly, varid_pdir_hourly, varid_pdif_hourly
    integer               ::  varid_skt_24h_history, varid_pdir_24h_history, varid_pdif_24h_history 
#endif
#ifdef with_online_nox
    integer               ::  varid_cp_daily, varid_lsp_daily
    integer               ::  varid_cp_history, varid_lsp_history
    integer               ::  varid_pulsing, varid_plsday, varid_plsdurat
#endif
#ifdef with_m7
    integer               ::  varid_h2o, varid_rw, varid_rwd !! , varid_ini_gph
    integer               ::  imode
    character(len=3), parameter :: h2o_name = 'h2o'
    character(len=3), parameter :: rwd_name = 'rwd'
    character(len=2), parameter :: rw_name  = 'rw'
#endif
    ! global work arrays to read data
    real, allocatable :: tmp4d(:,:,:,:)
    real, allocatable :: tmp3d(:,:,:) 
    real, allocatable :: rmt(:,:,:,:),rms(:,:,:,:), rmx(:,:,:,:),rmy(:,:,:,:), rmz(:,:,:,:)
#if defined(with_online_bvoc) || defined(with_online_nox)
    real, allocatable :: glb_sfc3D(:,:,:), glb_sfc4D(:,:,:,:)
    integer           :: n360, n180
#endif  

    ! for remapping:
    logical :: need_vremap, need_hremap, need_remap
    integer :: varid_at, varid_bt
    real    :: dx, dy
    real    :: factor_ch4
    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 ReadRc( rcF, 'restart.factor.ch4', factor_ch4, status, default = 1.0 )
    IF_ERROR_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( tmp4d(imr,jmr,lmr, ntracet) ) ! assume that ntracet is the max size needed
          allocate( tmp3d(imr,jmr,lmr+1       ) )
       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( tmp4d(1,1,1,1) )
          allocate( tmp3d(1,1,1)   )
       endif

       ! get extra work arrays for 1x1 dataset
#if defined(with_online_bvoc) || defined(with_online_nox)
       if(n==region1) then
          n360 = dgrid(iglbsfc)%im_region
          n180 = dgrid(iglbsfc)%jm_region
          if (isRoot) then
             allocate( glb_sfc3d(n360, n180, max(ndays_history, n_layers, ndrydays, nhours_history ) )
             allocate( glb_sfc4d(n360, n180, n_layers, ndays_history)    )
          else
             allocate( glb_sfc3d(1,1,1)   )
             allocate( glb_sfc4d(1,1,1,1) )
          endif
       end if
#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
          if ( do_m ) 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

#ifdef with_online_bvoc
          ! MEGAN/PCEEA history records; only once is ok
          if ( do_megan .and. (n==region1) ) then
             call MDF_Inq_VarID( ncid, 'skt_daily'       , varid_skt_daily       , status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Inq_VarID( ncid, 'skt_10d_history' , varid_skt_10d_history , status )
             IF_NOTOK_MDF(fid=ncid)
             if (megan) then
                call MDF_Inq_VarID( ncid, 'pdir_daily'       , varid_pdir_daily       , status )
                IF_NOTOK_MDF(fid=ncid)
                call MDF_Inq_VarID( ncid, 'pdir_10d_history' , varid_pdir_10d_history , status )
                IF_NOTOK_MDF(fid=ncid)
                call MDF_Inq_VarID( ncid, 'pdif_daily'       , varid_pdif_daily       , status )
                IF_NOTOK_MDF(fid=ncid)
                call MDF_Inq_VarID( ncid, 'pdif_10d_history' , varid_pdif_10d_history , status )
                IF_NOTOK_MDF(fid=ncid)
                call MDF_Inq_VarID( ncid, 'skt_hourly'       , varid_skt_hourly       , status )
                IF_NOTOK_MDF(fid=ncid)
                call MDF_Inq_VarID( ncid, 'skt_24h_history'  , varid_skt_24h_history  , status )
                IF_NOTOK_MDF(fid=ncid)
                call MDF_Inq_VarID( ncid, 'pdir_hourly'      , varid_pdir_hourly      , status )
                IF_NOTOK_MDF(fid=ncid)
                call MDF_Inq_VarID( ncid, 'pdir_24h_history' , varid_pdir_24h_history , status )
                IF_NOTOK_MDF(fid=ncid)
                call MDF_Inq_VarID( ncid, 'pdif_hourly'      , varid_pdif_hourly      , status )
                IF_NOTOK_MDF(fid=ncid)
                call MDF_Inq_VarID( ncid, 'pdif_24h_history' , varid_pdif_24h_history , status )
                IF_NOTOK_MDF(fid=ncid)
             else if (pceea) then
                call MDF_Inq_VarID( ncid, 'ssr_daily'        , varid_ssr_daily        , status )
                IF_NOTOK_MDF(fid=ncid)
                call MDF_Inq_VarID( ncid, 'ssr_10d_history'  , varid_ssr_10d_history  , status )
                IF_NOTOK_MDF(fid=ncid)
             endif
          end if
#endif
          
#ifdef with_online_nox
          ! precipitation history and pulsing parameters; only once is ok
          if ( do_pulse .and. (n==region1) ) then
             call MDF_Inq_VarID( ncid, 'cp_daily'   , varid_cp_daily,    status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Inq_VarID( ncid, 'lsp_daily'  , varid_lsp_daily,   status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Inq_VarID( ncid, 'cp_history' , varid_cp_history,  status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Inq_VarID( ncid, 'lsp_history', varid_lsp_history, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Inq_VarID( ncid, 'pulsing'    , varid_pulsing,     status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Inq_VarID( ncid, 'plsday'     , varid_plsday,      status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Inq_VarID( ncid, 'plsdurat'   , varid_plsdurat,    status )
             IF_NOTOK_MDF(fid=ncid)
          end if
#endif

#ifdef with_m7
#ifndef without_chemistry
          if (do_m) then
             ! M7 fields for optics
             call MDF_Inq_VarID( ncid, trim(h2o_name) , varid_h2o, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Inq_VarID( ncid, trim(rwd_name) , varid_rwd, status )
             IF_NOTOK_MDF(fid=ncid)
             call MDF_Inq_VarID( ncid, trim(rw_name ) , varid_rw, status  )
             IF_NOTOK_MDF(fid=ncid)
          end if
#endif
#endif
       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

       if ( do_m ) then
          write (gol,'("    restore air mass ...")'); call goPr
          
          if (isRoot) call MDF_Get_Var( ncid, varid_m, tmp3d(:,:,1:lmr), status)
          IF_NOTOK_MDF(fid=ncid)
          
          call scatter( dgrid(n), m_dat(n)%data, tmp3d(:,:,1:lmr), m_dat(n)%halo, status)
          IF_NOTOK_RETURN(status=1)
       end if

       !! surface fluxes
       !if ( do_sflux ) then
       !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.e30")') 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.e30")') 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)
                         call Fill3D( global_lli(n), levi, 'n', sp_gbl(:,:,1), rmt(:,:,:,itr), &
                                      lli_restart,   levi_restart, src_glb, 'sum', status )
                         IF_NOTOK_RETURN(status=1)
                      else
                         call MDF_Get_Var( ncid, varid_rm, rmt(:,:,:,itr), status, start=(/1,1,1,itr_file/))
                         IF_NOTOK_MDF(fid=ncid)
                      endif

                      ! Scale methane concentrations by a factor specified in the rc file
                      if ( (factor_ch4 /= 1.) .and. (itr == ich4) ) then
                         write(gol,*) '*WARNING*: CH4 mixing ratio and slopes from restart file'; call goPr
                         write(gol,*) '*WARNING*: ... scaled by a factor: ', factor_ch4; call goPr
                         rmt(:,:,:,itr) = rmt(:,:,:,itr) * factor_ch4
                      endif
#ifdef slopes
                      ! read slopes
                      if (.not. need_remap) 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)
                         
                         ! Scale methane concentration slopes by a factor specified in the rc file
                         if ( (factor_ch4 /= 1.) .and. (itr == ich4) ) then
                            mass_dat(n)%rxm(:,:,:,itr)= mass_dat(n)%rxm(:,:,:,itr) * factor_ch4
                            mass_dat(n)%rym(:,:,:,itr)= mass_dat(n)%rym(:,:,:,itr) * factor_ch4
                            mass_dat(n)%rzm(:,:,:,itr)= mass_dat(n)%rzm(:,:,:,itr) * factor_ch4
                         endif
                      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)
                         call Fill3D( global_lli(n), levi, 'n', sp_gbl(:,:,1), rms(:,:,:,itr), &
                                      lli_restart,   levi_restart, src_glb, 'sum', status )
                         IF_NOTOK_RETURN(status=1)
                      else
                         call MDF_Get_Var( ncid, varid_rmc, rms(:,:,:,itr), status, start=(/1,1,1,itr_file/) )
                         IF_NOTOK_MDF(fid=ncid)
                      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) 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 "READRM"
       deallocate(rmt)
       if ( ntrace_chem > 0 ) deallocate(rms)
#ifdef slopes
       deallocate(rmx, rmy, rmz)
#endif       

#ifdef with_online_bvoc
       ! MEGAN/PCEEA history records; read only once
       if ( do_megan .and. (n==region1) ) then

          if (isRoot) call MDF_Get_Var( ncid, varid_skt_daily, glb_sfc3d(:,:,1), status)
          IF_NOTOK_MDF(fid=ncid)
          call scatter( dgrid(n), skt_daily, glb_sfc3d(:,:,1), 0, status)

          
          if (isRoot) call MDF_Get_Var( ncid, varid_skt_10d_history, glb_sfc3d(:,:, 1:ndays_history), status)
          IF_NOTOK_MDF(fid=ncid)
          call scatter( dgrid(n), skt_10d_history(:,:, 1:ndays_history), glb_sfc3d(:,:, 1:ndays_history), 0, status)

          if (megan) then
             write (gol,'("  restore MEGAN history parameters ...")'); call goPr

             if (isRoot) call MDF_Get_Var( ncid, varid_pdir_daily, glb_sfc3d(:,:,1), status)
             IF_NOTOK_MDF(fid=ncid)
             call scatter( dgrid(n), pdir_daily, glb_sfc3d(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)
             
             if (isRoot) call MDF_Get_Var( ncid, varid_pdir_10d_history, glb_sfc3d(:,:, 1:ndays_history), status)
             IF_NOTOK_MDF(fid=ncid)
             call scatter( dgrid(n), pdir_10d_history(:,:, 1:ndays_history), glb_sfc3d(:,:, 1:ndays_history), 0, status)
             IF_NOTOK_RETURN(status=1)
             
             if (isRoot) call MDF_Get_Var( ncid, varid_pdif_daily, glb_sfc3d (:,:,1:n_layers), status)
             IF_NOTOK_MDF(fid=ncid)
             call scatter( dgrid(n), pdif_daily, glb_sfc3d(:,:,1:n_layers), 0, status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) call MDF_Get_Var( ncid, varid_pdif_10d_history, glb_sfc4D, status)
             IF_NOTOK_MDF(fid=ncid) 
             call scatter( dgrid(n), pdif_10d_history, glb_sfc4D, 0, status)
             IF_NOTOK_RETURN(status=1)
             
             if (isRoot) call MDF_Get_Var( ncid, varid_skt_hourly, glb_sfc3d(:,:,1), status)
             IF_NOTOK_MDF(fid=ncid)
             call scatter( dgrid(n), skt_hourly, glb_sfc3d(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)
             
             if (isRoot) call MDF_Get_Var( ncid, varid_skt_24h_history, glb_sfc3d(:,:, 1:nhours_history), status)
             IF_NOTOK_MDF(fid=ncid)
             call scatter( dgrid(n), skt_24h_history, glb_sfc3d(:,:,1:nhours_history), 0, status)
             IF_NOTOK_RETURN(status=1)
             
             if (isRoot) call MDF_Get_Var( ncid, varid_pdir_hourly, glb_sfc3d(:,:,1), status)
             IF_NOTOK_MDF(fid=ncid)
             call scatter( dgrid(n), pdir_hourly, glb_sfc3d(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)
             
             if (isRoot) call MDF_Get_Var( ncid, varid_pdir_24h_history, glb_sfc3d(:,:, 1:nhours_history), status)
             IF_NOTOK_MDF(fid=ncid)
             call scatter( dgrid(n), pdir_24h_history, glb_sfc3d(:,:,1:nhours_history), 0, status)
             IF_NOTOK_RETURN(status=1)
             
             if (isRoot) call MDF_Get_Var( ncid, varid_pdif_hourly, glb_sfc3d(:,:, 1:n_layers), status)
             IF_NOTOK_MDF(fid=ncid)
             call scatter( dgrid(n), pdif_hourly, glb_sfc3d(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)
             
             if (isRoot) call MDF_Get_Var( ncid, varid_pdif_24h_history, glb_sfc4D, status)
             IF_NOTOK_MDF(fid=ncid)
             call scatter( dgrid(n), pdif_24h_history, glb_sfc4D, 0, status)
             IF_NOTOK_RETURN(status=1)

          else if (pceea) then
             write (gol,'("  restore PCEEA history parameters ...")'); call goPr

             if (isRoot) call MDF_Get_Var( ncid, varid_ssr_daily, glb_sfc3d(:,:,1), status)
             IF_NOTOK_MDF(fid=ncid)
             call scatter( dgrid(n), ssr_daily, glb_sfc3d(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)
             
             if (isRoot) call MDF_Get_Var( ncid, varid_ssr_10d_history,  ssr_10d_history(:,:, 1:ndays_history), status)
             IF_NOTOK_MDF(fid=ncid)
             call scatter( dgrid(n), ssr_10d_history(:,:, 1:ndays_history), glb_sfc3d(:,:, 1:ndays_history), 0, status)
             IF_NOTOK_RETURN(status=1)

          endif
       end if
#endif

#ifdef with_online_nox
       ! precipitation history and pulsing parameters; read only once
       if ( do_pulse .and. (n==region1) ) then
          
          write (gol,'("  restore precipitation history and pulsing parameters ...")'); call goPr

          if (isRoot) call MDF_Get_Var( ncid,  varid_cp_daily, glb_sfc3d(:,:,1), status)
          IF_NOTOK_MDF(fid=ncid)
          call scatter( dgrid(n), cp_daily, glb_sfc3d(:,:,1), 0, status)
          IF_NOTOK_RETURN(status=1)
          
          if (isRoot) call MDF_Get_Var( ncid,  varid_lsp_daily, glb_sfc3d(:,:,1), status)
          IF_NOTOK_MDF(fid=ncid)
          call scatter( dgrid(n), lsp_daily, glb_sfc3d(:,:,1), 0, status)
          IF_NOTOK_RETURN(status=1)
          
          if (isRoot) call MDF_Get_Var( ncid,  varid_cp_history,  glb_sfc3d(:,:, 1:ndrydays), status)
          IF_NOTOK_MDF(fid=ncid)
          call scatter( dgrid(n), cp_history(:,:, 1:ndrydays), glb_sfc3d(:,:, 1:ndrydays), 0, status)
          IF_NOTOK_RETURN(status=1)
          
          if (isRoot) call MDF_Get_Var( ncid,  varid_lsp_history, glb_sfc3d(:,:, 1:ndrydays), status)
          IF_NOTOK_MDF(fid=ncid)
          call scatter( dgrid(n), lsp_history(:,:, 1:ndrydays), glb_sfc3d(:,:, 1:ndrydays), 0, status)
          IF_NOTOK_RETURN(status=1)

          if (isRoot) call MDF_Get_Var( ncid,  varid_pulsing, glb_sfc3d(:,:,1), status)
          IF_NOTOK_MDF(fid=ncid)
          call scatter( dgrid(n), pulsing_field, glb_sfc3d(:,:,1), 0, status)
          IF_NOTOK_RETURN(status=1)

          if (isRoot) call MDF_Get_Var( ncid,  varid_plsday, glb_sfc3d(:,:,1), status)
          IF_NOTOK_MDF(fid=ncid)
          call scatter( dgrid(n), plsday_field, glb_sfc3d(:,:,1), 0, status)
          IF_NOTOK_RETURN(status=1)

          if (isRoot) call MDF_Get_Var( ncid,  varid_plsdurat, glb_sfc3d(:,:,1), status)
          IF_NOTOK_MDF(fid=ncid)
          call scatter( dgrid(n), plsdurat_field, glb_sfc3d(:,:,1), 0, status)
          IF_NOTOK_RETURN(status=1)

       end if
#endif

#ifdef with_m7
#ifndef without_chemistry
       
       if (do_m) then

          write (gol,'("    restore M7 fields for optics ...")'); call goPr

          ! water: get 4d array
          if (isRoot) call MDF_Get_Var( ncid, varid_h2o, tmp4d(:,:,:,1:nsol), status )
          IF_NOTOK_MDF(fid=ncid)

          do imode=1,nsol
             call scatter( dgrid(n), h2o_mode(n,imode)%d3, tmp4d(:,:,:,imode), h2o_mode(n,imode)%halo, status)
          end do

          ! dry radii: get 4d array 
          if (isRoot) call MDF_Get_Var( ncid, varid_rwd, tmp4d(:,:,:,1:nsol), status)
          IF_NOTOK_MDF(fid=ncid)
       
          do imode=1,nsol
             call scatter( dgrid(n), rwd_mode(n,imode)%d3, tmp4d(:,:,:,imode), rwd_mode(n,imode)%halo, status)
          end do
          
          ! (wet) radii: get 4d array 
          if (isRoot) call MDF_Get_Var( ncid, varid_rw, tmp4d(:,:,:,1:nmod), status)
          IF_NOTOK_MDF(fid=ncid)
          
          do imode=1,nmod
             call scatter( dgrid(n), rw_mode(n,imode)%d3, tmp4d(:,:,:,imode), rw_mode(n,imode)%halo, status)
          end do
       end if
#endif
#endif /* with_m7 */

       if (isRoot) call MDF_Close( ncid, status )
       IF_NOTOK_RETURN(status=1)
       
       deallocate( tmp4d )
       deallocate( tmp3d )
#if defined(with_online_bvoc) || defined(with_online_nox)
       deallocate(glb_sfc3D, glb_sfc4D)
#endif  

    ENDDO REG

    status = 0

  END SUBROUTINE RESTART_READ
  !EOC

END MODULE RESTART