ecearth3/sources/tm5mp/base/io_save.F90

!### macro's #####################################################
!
#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_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
#define IF_NOTOK_MDF(action) if (status/=0) then; TRACEBACK; action; call MDF_CLose(fid,status); status=1; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5-MP                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      IO_SAVE
!
! !DESCRIPTION: Contains routines to read HDF emissions files and 
!               the main model state (as replacement for restart file).
!\\
!\\
! !INTERFACE: 
!
MODULE IO_SAVE
  !
  ! !USES:
  !
  use GO,          only : gol, goErr, goPr, goBug
  use GO,          only : TrcFile, InitRc => Init, DoneRc => Done, ReadRc
  use global_data, only : rcfile

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  PUBLIC :: READHDFMMR        ! Read mmix from "old save file"                                      (istart=4 )
  PUBLIC :: READ_MMIX         ! Read  mmix file (mean mixing ratio of transported tracers)          (istart=5 )
  PUBLIC :: READ_SAVE_FILE_30 ! Read  save-file (mixing ratio, opt. slopes, of transported tracers) (istart=30)
  PUBLIC :: READ_SAVE_FILE    ! Read  save file (mass of all tracers, grid info)                    (istart=31)
  PUBLIC :: WRITE_SAVE_FILE   ! Write save-file (to be read with read_save_file)

  PUBLIC :: READ_HDF4         ! Read a record from hdf file (new way)
  PUBLIC :: READTM3HDF        ! Read a record from hdf file (old way, a la TM3)
  PUBLIC :: READ_SCATTER_HDF  ! Read a record from hdf file, and do the MPI scattering
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  :: mname = 'io_save'
  type(TrcFile)                :: rcF
  integer                      :: fid  ! file id for IF_NOTOK_MDF macro
  !
  ! !INTERFACE:
  !
  INTERFACE READ_HDF4
     MODULE PROCEDURE READ_HDF4_2DR
  END INTERFACE
  !
  ! !REVISION HISTORY: 
  !   16 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !  (1) The READTM3 routines are kept for quick fix when converting old code.
  !      All new code should use the MDF module to read any data. 
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    READ_HDF4_2DR
  !
  ! !DESCRIPTION: read one 2D record of real from an HDF4 file.
  ! 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE READ_HDF4_2DR(fname, odata, status, dir, vname, vid)
    !
    ! !USES:
    !
    USE MDF,         ONLY : MDF_Open, MDF_HDF4, MDF_READ, MDF_Inq_VarID, MDF_Get_Var, MDF_Close
    USE global_data, ONLY : inputdir
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*), intent(in)           :: fname  ! name of HDF file to read
    character(len=*), intent(in), optional :: dir    ! dir where file is. Default is "inputdir"
    character(len=*), intent(in), optional :: vname  ! record name. Default is to use record #0.
    integer,          intent(in), optional :: vid    ! record id #. Default is to use record #0.
    !
    ! !OUTPUT PARAMETERS:
    !
    real,            intent(out) :: odata(:,:)
    integer,         intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !   28 Mar 2012 - P. Le Sager - v0
    !
    ! !REMARKS:
    !  (1) This wrapper around MDF routines is only for convenience: a lot of
    !       calls for the chemistry emissions, and participate to the effort
    !       of getting rid of file_hdf.f90 routine, in favor of the MDF module.
    !  (2) By default **FIRST** record is read, unless another variable number
    !       (vid) and/or variable name (vname) is passed. If both vid and vname
    !       are passed, vname is ignored.
    !  
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/read_hdf4_2dr'
    character(len=512):: Lvname, Ldir
    integer           :: Lvid
 
    !--- check input
    
    Ldir=inputdir
    if(present(dir)) Ldir=dir

    Lvname=''
    if(present(vname)) Lvname=vname

    Lvid=1
    if (present(vid)) then
       Lvid=vid
       Lvname=''
    end if

    !--- read
    
    CALL MDF_Open( trim(Ldir)//'/'//TRIM(fname), MDF_HDF4, MDF_READ, fid, status )
    IF_NOTOK_RETURN(status=1)

    if (len(trim(Lvname))/=0) then
       CALL MDF_Inq_VarID( fid, TRIM(Lvname), Lvid, status )
       IF_NOTOK_MDF(status=1)
    end if
    
    CALL MDF_Get_Var( fid, Lvid, odata, status )
    IF_NOTOK_MDF(status=1)

    CALL MDF_Close( fid, status )
    IF_NOTOK_RETURN(status=1)

    status=0

  END SUBROUTINE READ_HDF4_2DR
  !EOC

  ! ======================================================================


  !WP! this subroutine reads HDF emissions files on 
  !    1x1 degree resolution and returns the specified field.
  !WP! input for this subroutine are:
  !WP!
  !WP! filename --> name of hdf file to be opened or read
  !MK! rank    --> rank of the dataset (1,2,3)
  !WP! im_emis --> x-DIMENSION of requested field
  !WP! jm_emis --> y-DIMENSION of requested field     
  !WP! lm_emis --> z-DIMENSION of requested field
  !WP! month  --> the month to read (0 for january or simply the first field)
  !WP! field --> the field (im_emis,jm_emis,lm_emis) to put the values in
  !WP! field_name --> the name of the field to match the one in the HDF file
  !WP!
  !WP! The routine checks the dimensions and reports the succes/failure.
  
  subroutine readtm3hdf( filename, rank, im_emis, jm_emis, lm_emis, &
       month, field, field_name, idt )

    use file_hdf, only : THdfFile, Init, Done
    use io_hdf  , only : FAIL
    use io_hdf  , only : io_read2d_32g, io_read3d_32g
    use toolbox , only : escape_tm

    ! --- in/out --------------------------------------

    character(len=*),intent(in)  :: filename
    character(len=*),intent(in)  :: field_name
    integer,intent(in)           :: rank
    integer,intent(in)           :: month
    integer,intent(in)           :: im_emis, jm_emis, lm_emis
    real,dimension(im_emis,jm_emis,lm_emis),intent(out)  :: field
    integer,dimension(6),optional :: idt

    ! --- local ----------------------------------------

    type(THdfFile)           ::  hdf
    integer                  ::  status
    integer                  ::  idx

    real, allocatable        ::  field1d(:)   ! first dim will be 1
    real, allocatable        ::  field2d(:,:)
    real, allocatable        ::  field3d(:,:,:)

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

    ! open hdf file:    
    call Init( hdf, trim(filename), 'read', status )
    if (status/=0) call escape_tm('ERROR in readtm3hdf')

    !FD hdf files indices run from 0 ...n-1
    idx = month

    select case(rank) 
    case(1)    ! latitudional field

       allocate(field1d(jm_emis))
       if(present(idt)) call io_read2d_32g(hdf%id,im_emis,jm_emis,field1d,field_name, status,idate=idt)
       if(.not.present(idt)) call io_read2d_32g(hdf%id,im_emis,jm_emis,field1d,field_name, status,index=idx)
       if ( status /= 0 ) then
          write(gol,*) 'readtm3hdf: io_read2d_32g failed:'            ; call goPr
          write(gol,*) 'readtm3hdf:   filename   : ', trim(filename)  ; call goPr
          write(gol,*) 'readtm3hdf:   rank       : ', rank            ; call goPr
          write(gol,*) 'readtm3hdf:   im_emis    : ', im_emis         ; call goPr
          write(gol,*) 'readtm3hdf:   jm_emis    : ', jm_emis         ; call goPr
          write(gol,*) 'readtm3hdf:   lm_emis    : ', lm_emis         ; call goPr
          write(gol,*) 'readtm3hdf:   month      : ', month           ; call goPr
          write(gol,*) 'readtm3hdf:   field_name : ', trim(field_name); call goPr
          if (present(idt)) write(gol,*) 'readtm3hdf:   idt        : ', idt; call goPr
          call escape_tm('readtm3hdf: Error reading HDF file 1D')
       end if
       field(1,:,1)=field1d
       deallocate(field1d)

    case(2)    ! 2D surface field

       allocate(field2d(im_emis,jm_emis))
       if(present(idt)) call io_read2d_32g(hdf%id,im_emis,jm_emis,field2d,field_name, status,idate=idt)
       if(.not.present(idt)) call io_read2d_32g(hdf%id,im_emis,jm_emis,field2d,field_name, status,index=idx)
       if ( status /= 0 ) then
          write(gol,*) 'readtm3hdf: io_read2d_32g failed:'                        ; call goPr
          write(gol,*) 'readtm3hdf:   filename   : ', trim(filename)              ; call goPr
          write(gol,*) 'readtm3hdf:   rank       : ', rank                        ; call goPr
          write(gol,*) 'readtm3hdf:   im_emis    : ', im_emis                     ; call goPr
          write(gol,*) 'readtm3hdf:   jm_emis    : ', jm_emis                     ; call goPr
          write(gol,*) 'readtm3hdf:   lm_emis    : ', lm_emis                     ; call goPr
          write(gol,*) 'readtm3hdf:   month      : ', month                       ; call goPr
          write(gol,*) 'readtm3hdf:   field_name : ', trim(field_name)            ; call goPr
          if (present(idt)) write(gol,*) 'readtm3hdf:   idt        : ', idt       ; call goPr
          call escape_tm('readtm3hdf: Error reading HDF file 2D')
       endif
       field(:,:,1)=field2d
       deallocate(field2d)

    case (3)   ! 3D field

       allocate(field3d(im_emis,jm_emis,lm_emis))
       if(present(idt)) call io_read3d_32g(hdf%id,im_emis,jm_emis,lm_emis,field3d,field_name, status,idate=idt)
       if(.not.present(idt)) call io_read3d_32g(hdf%id,im_emis,jm_emis,lm_emis,field3d,field_name, status,index=idx)
       if ( status /= 0 )  then
          write(gol,*) 'readtm3hdf: io_read2d_32g failed:'                   ; call goPr
          write(gol,*) 'readtm3hdf:   filename   : ', trim(filename)         ; call goPr
          write(gol,*) 'readtm3hdf:   rank       : ', rank                   ; call goPr
          write(gol,*) 'readtm3hdf:   im_emis    : ', im_emis                ; call goPr
          write(gol,*) 'readtm3hdf:   jm_emis    : ', jm_emis                ; call goPr
          write(gol,*) 'readtm3hdf:   lm_emis    : ', lm_emis                ; call goPr
          write(gol,*) 'readtm3hdf:   month      : ', month                  ; call goPr
          write(gol,*) 'readtm3hdf:   field_name : ', trim(field_name)       ; call goPr
          if (present(idt)) write(gol,*) 'readtm3hdf:   idt        : ', idt  ; call goPr
          call escape_tm('readtm3hdf: Error reading HDF file 3D')
       endif
       field(:,:,:)=field3d
       deallocate(field3d)

    case (23)   ! 2D field lat-pres

       allocate(field2d(jm_emis,lm_emis))
       if(present(idt)) call io_read2d_32g(hdf%id,jm_emis,lm_emis,field2d,field_name, status,idate=idt)
       if(.not.present(idt)) call io_read2d_32g(hdf%id,jm_emis,lm_emis,field2d,field_name, status,index=idx)
       if ( status /= 0 )  then   !might be due to (1,jm,lm) file!
          write(gol,*) 'readtm3hdf: try shape : ', 1, jm_emis, lm_emis
          allocate(field3d(1,jm_emis,lm_emis))
          if ( present(idt) ) then
             call io_read3d_32g(hdf%id,1,jm_emis,lm_emis,field3d,field_name, status,idate=idt)
          else
             call io_read3d_32g(hdf%id,1,jm_emis,lm_emis,field3d,field_name, status,index=idx)
          end if
          if ( status /= 0 ) then
             write(gol,*) 'readtm3hdf: io_read2d_32g failed:'                   ; call goPr
             write(gol,*) 'readtm3hdf:   filename   : ', trim(filename)         ; call goPr
             write(gol,*) 'readtm3hdf:   rank       : ', rank                   ; call goPr
             write(gol,*) 'readtm3hdf:   im_emis    : ', im_emis                ; call goPr
             write(gol,*) 'readtm3hdf:   jm_emis    : ', jm_emis                ; call goPr
             write(gol,*) 'readtm3hdf:   lm_emis    : ', lm_emis                ; call goPr
             write(gol,*) 'readtm3hdf:   month      : ', month                  ; call goPr
             write(gol,*) 'readtm3hdf:   field_name : ', trim(field_name)       ; call goPr
             if (present(idt)) write(gol,*) 'readtm3hdf:   idt        : ', idt  ; call goPr
             call escape_tm('readtm3hdf: Error reading HDF file 23')
          end if
          field(:,:,:)=field3d
          deallocate(field3d)
       else
          field(1,:,:)=field2d
       end if
       deallocate(field2d)

    case default
       write (gol,*) 'readtm3hdf: Please check the rank of the required field,', &
            ' must be 1,2,3, or 23'; call goPr
    end select

    ! close:
    call Done( hdf, status )
    if (status/=0) call escape_tm('ERROR in readtm3hdf')

    !if ( idx >= 0 ) then
    !  write (gol,*) 'readtm3hdf: ',trim(filename),' ; ds#= ',idx+1; call goPr
    !  write (gol,*) 'readtm3hdf: total month sum after reading the full field', &
    !     sum(field)/1e9 ,' (x1e9) '; call goPr
    !end if

  end subroutine readtm3hdf
  

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    READ_SCATTER_HDF
  !
  ! !DESCRIPTION: simple adaptation of readtm3hdf to use MDF module, and
  !               scatter data.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE READ_SCATTER_HDF( filename, rank, region, i1, j1, field, status, &
                               idx,  field_name, idt )
    !
    ! !USES:
    !
    use partools,      only : isRoot
    use tm5_distgrid,  only : dgrid, Get_DistGrid, scatter
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*), intent(in)  :: filename            ! HDF4 file to read
    integer,          intent(in)  :: rank                ! 1,2,.. [only 2 implemented for now]
    integer,          intent(in)  :: region, i1, j1      ! region and local array start index

    integer,          intent(in), optional :: idx        ! variable id in [1..nmax]
    character(len=*), intent(in), optional :: field_name ! variable name, ignored if idx set 
    integer,          intent(in), optional :: idt(6)     ! date [not implemented]
    !
    ! !OUTPUT PARAMETERS:
    !
    real,    intent(out)  :: field(i1:,j1:)  ! scattered read data
    integer, intent(out)  :: status
    !
    ! !REVISION HISTORY: 
    !   23 Jan 2012 - P. Le Sager - for lat-lon mpi decomposition
    !
    ! !REMARKS:
    !   (1) tested only with rank=2 (see dry_deposition.F90)
    !   (2) not all features of readtm3hdf are available, but probably also
    !        not needed... will add as needed.
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer           :: im_emis, jm_emis                      
    real, allocatable :: field1d(:)
    real, allocatable :: field2d(:,:)

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

    CALL GET_DISTGRID( DGRID(region), I_WRLD=im_emis, J_WRLD=jm_emis )

    status=0
    
    READ: if (isRoot) then

       allocate(field2d(im_emis,jm_emis))
       field2D=0.

       select case(rank)

!       case(1)    ! latitudional field
!
!           allocate(field1d(jm_emis))
!           if(present(idt)) call io_read2d_32g(hdf%id,im_emis,jm_emis,field1d,field_name, status,idate=idt)
!           if(.not.present(idt)) call io_read2d_32g(hdf%id,im_emis,jm_emis,field1d,field_name, status,index=idx)
!           if ( status /= 0 ) then
!              write(gol,*) 'readtm3hdf: io_read2d_32g failed:'
!              write(gol,*) 'readtm3hdf:   filename   : ', trim(filename)
!              write(gol,*) 'readtm3hdf:   rank       : ', rank
!              write(gol,*) 'readtm3hdf:   im_emis    : ', im_emis
!              write(gol,*) 'readtm3hdf:   jm_emis    : ', jm_emis
!              write(gol,*) 'readtm3hdf:   month      : ', month
!              write(gol,*) 'readtm3hdf:   field_name : ', trim(field_name)
!              if (present(idt)) write(gol,*) 'readtm3hdf:   idt        : ', idt
!              call escape_tm('readtm3hdf: Error reading HDF file 1D')
!           end if
!           field2D(1,:)=field1d
!           deallocate(field1d)

       case(2)    ! 2D surface field


          call READ_HDF4(filename, field2d, status, dir='', vname=field_name, vid=idx)
          IF_NOTOK_RETURN(status=1)
          

       case default
          write (gol,*) rname//' : Please check the rank of the required field,', &
               ' must be 1,2,3, or 23'; call goPr
          status=1
       end select

    else

       allocate(field2d(1,1))

    end if READ

    call SCATTER( dgrid(region), field, field2d, 0, status)
    IF_NOTOK_RETURN(status=1)

    deallocate(field2d)

  END SUBROUTINE READ_SCATTER_HDF
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    READ_SAVE_FILE_30
  !
  ! !DESCRIPTION: reads initial MIXING RATIO (+ slopes, and second moments if
  !               they are switched on) of TRANSPORTED tracers from HDF-4
  !               "save files". The save file must be specified in rc file with key:
  !               
  !               start.30.<region-name>
  !               
  !               No regridding available.
  !
  !               Used for istart=30.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE READ_SAVE_FILE_30( region, file_name, status )
    !
    ! !USES:
    !
    use io_hdf
    use dims,          only : im, jm, lm
    use global_data,   only : mass_dat
    use meteodata  ,   only : m_dat
    use chem_param,    only : ntracet
    use partools,      only : isRoot
    use tm5_distgrid,  only : dgrid, Get_DistGrid, scatter
    !
    ! !INPUT PARAMETERS:
    !
    integer,          intent(in) :: region
    character(len=*), intent(in) :: file_name 
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,         intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !    4 Oct 2011 - P. Le Sager - TM5-MP version
    !
    ! !REMARKS:
    ! (1) no zoom (must put from_file stuff back in)
    ! (2) FIXME: not tested...
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    real,dimension(:,:,:,:),pointer :: rm
    real,dimension(:,:,:,:),pointer :: rxm,rym,rzm
#ifdef secmom
    real,dimension(:,:,:,:),pointer :: rxxm,rxym,rxzm,ryym,ryzm,rzzm
#endif
    real,dimension(:,:,:),pointer   :: m

    integer                         :: io, sfstart
    integer                         :: istat
    integer                         :: sfend
    integer                         :: i,j,l,n
    integer                         :: imr,jmr,lmr, halo
    real,dimension(:,:,:,:), allocatable :: field
    real,dimension(:,:,:),   allocatable :: msave, msave_local

    integer :: i1,  i2,  j1,  j2

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

    m  => m_dat(region)%data
    rm => mass_dat(region)%rm

#ifdef slopes
    rxm => mass_dat(region)%rxm
    rym => mass_dat(region)%rym
    rzm => mass_dat(region)%rzm
#ifdef secmom
    rxxm => mass_dat(region)%rxxm
    rxym => mass_dat(region)%rxym
    rxzm => mass_dat(region)%rxzm
    ryym => mass_dat(region)%ryym
    ryzm => mass_dat(region)%ryzm
    rzzm => mass_dat(region)%rzzm
#endif
#endif

    halo = mass_dat(region)%halo
    
    ! global array
    imr = im(region) ; jmr = jm(region) ; lmr = lm(region)

    if ( isRoot ) then
       allocate(field(imr, jmr, lmr, ntracet))
       field = 0.0
       allocate(msave(imr,jmr,lmr))
    else
       allocate(field(1,1,1,1))
       allocate(msave(1,1,1))
    end if


    ! local indices
    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2  )
    allocate(msave_local(i1:i2,j1:j2,lmr))
    
    
    ! read on root
    if ( isRoot ) then

       io = sfstart( file_name, DFACC_READ )

       if ( io == FAIL ) then 
          write(gol,*) 'Read_save_file_30 : Could not open file '// trim(file_name) ; call goPr
          call goErr
          TRACEBACK ; status=1 ; return
       end if

       write(gol,'("Read_save_file_30: opened for READ access: ",a)') trim(file_name); call goPr

       call io_read3d_32g(io,imr,jmr,lmr,msave,'m',status)
       IF_NOTOK_RETURN(status=1)

       call io_read4d_32g(io,imr,jmr,lmr,ntracet, field,'rm',status)

       ! try again
       if( status /= 0 .and. ntracet == 1 ) then
          call io_read3d_32g(io,imr,jmr,lmr,field,'rm',status)
       end if
       IF_NOTOK_RETURN(write (*,'("read_save_file_30: Read rm, status = ",i6)') status)
       
    end if  !root


    ! get local chunk
    call SCATTER( dgrid(region), rm, field, halo, status)
    IF_NOTOK_RETURN(status=1)
    
    field = 0.0

    call SCATTER( dgrid(region), msave_local, msave, 0, status)
    IF_NOTOK_RETURN(status=1)
    
#ifdef slopes

    if ( isRoot ) then

       call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rxm',status)

       if ( status /= 0 .and. ntracet == 1 ) then
          call io_read3d_32g(io,imr,jmr,lmr,field,'rxm',status)
       end if
       IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status)

    end if

    ! get local chunk
    call SCATTER( dgrid(region), rxm, field, halo, status)
    IF_NOTOK_RETURN(status=1)

    field = 0.0


    if(isRoot) then

       call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rym',status)

       if ( status /= 0 .and. ntracet == 1 ) then
          call io_read3d_32g(io,imr,jmr,lmr,field,'rym',status)
       end if
       IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status)

    end if  !root

    ! get local chunk
    call SCATTER( dgrid(region), rym, field, halo, status)
    IF_NOTOK_RETURN(status=1)

    field = 0.0


    if ( isRoot ) then

       call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rzm',status)
       if ( status /= 0 .and. ntracet == 1 ) then
          call io_read3d_32g(io,imr,jmr,lmr,field,'rzm',status)
       end if
       IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status)

    end if

    ! get local chunk
    call SCATTER( dgrid(region), rzm, field, halo, status)
    IF_NOTOK_RETURN(status=1)

    ! reset
    field = 0.0

#endif


    ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    ! >>> second moments
    ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#ifdef secmom

    if ( isRoot) then

       call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rxxm',status)
       if ( status /= 0 .and. ntracet == 1 ) then
          call io_read3d_32g(io,imr,jmr,lmr,field,'rxxm',status)
       end if
       IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status)

    endif

    call SCATTER( dgrid(region), rxxm, field, halo, status)
    IF_NOTOK_RETURN(status=1)

    field = 0.0


    if(isRoot) then

       call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rxym',status)
       if ( status /= 0 .and. ntracet == 1 ) then
          call io_read3d_32g(io,imr,jmr,lmr,field,'rxym',status)
       end if
       IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status)

    endif

    call SCATTER( dgrid(region), rxym, field, halo, status)
    IF_NOTOK_RETURN(status=1)

    field = 0.0


    if(isRoot) then

       call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rxzm',status)
       if ( status /= 0 .and. ntracet == 1 ) then
          call io_read3d_32g(io,imr,jmr,lmr,field,'rxzm',status)
       end if
       IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status)

    endif

    call SCATTER( dgrid(region), rxzm, field, halo, status)
    IF_NOTOK_RETURN(status=1)

    field = 0.0


    if(isRoot) then

       call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'ryym',status)
       if ( status /= 0 .and. ntracet == 1 ) then
          call io_read3d_32g(io,imr,jmr,lmr,field,'ryym',status)
       end if
       IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status)

    endif

    call SCATTER( dgrid(region), ryym, field, halo, status)
    IF_NOTOK_RETURN(status=1)

    field = 0.0


    if(isRoot) then

       call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'ryzm',status)
       if ( status /= 0 .and. ntracet == 1 ) then
          call io_read3d_32g(io,imr,jmr,lmr,field,'ryzm',status)
       end if
       IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status)


    endif

    call SCATTER( dgrid(region), ryzm, field, halo, status)
    IF_NOTOK_RETURN(status=1)

    ! reset
    field = 0.0


    if(isRoot) then

       call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rzzm',status)
       if ( status /= 0 .and. ntracet == 1 ) then
          call io_read3d_32g(io,imr,jmr,lmr,field,'rzzm',status)
       end if
       IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status)

    endif

    call SCATTER( dgrid(region), rzzm, field, halo, status)
    IF_NOTOK_RETURN(status=1)

    field = 0.0
#endif
    

    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    ! <<< end second moments
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

    if ( isRoot ) then

       istat = sfend(io)
       if ( istat /= FAIL ) then
          write(gol,'("read_save_file_30: ... file closed")') ; call goPr
          write(gol,'("read_save_file_30: ")') ; call goPr
       else
          TRACEBACK ; status=1 ; return
       end if

    end if

    ! CONVERSION
    ! rm/msave ---> mixing ratio on file. mr*m --> something back in kilo's

    if ( isRoot ) then

       write(gol,*) 'read_save_file_30: Maximum ratio of the air masses:', &
            maxval(abs(msave/m(1:imr,1:jmr,1:lmr))) ; call goPr
       write(gol,*) 'read_save_file_30: Minimum ratio of the air masses:', &
            minval(abs(msave/m(1:imr,1:jmr,1:lmr))) ; call goPr

    end if

    do n = 1, ntracet
       do l = 1, lmr 
          rm(i1:i2,j1:j2,l,n)  = rm (i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l)
#ifdef slopes
          rxm(i1:i2,j1:j2,l,n) = rxm(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l)
          rym(i1:i2,j1:j2,l,n) = rym(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l)
          rzm(i1:i2,j1:j2,l,n) = rzm(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l)
#ifdef secmom
          rxxm(i1:i2,j1:j2,l,n) = rxxm(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l)
          rxym(i1:i2,j1:j2,l,n) = rxym(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l)
          rxzm(i1:i2,j1:j2,l,n) = rxzm(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l)
          ryym(i1:i2,j1:j2,l,n) = ryym(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l)
          ryzm(i1:i2,j1:j2,l,n) = ryzm(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l)
          rzzm(i1:i2,j1:j2,l,n) = rzzm(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l)
#endif
#endif
       end do
    end do

    ! cleanup
    deallocate(msave, msave_local)
    deallocate(field)

    nullify(m)
    nullify(rm)
#ifdef slopes
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
#ifdef secmom
    nullify(rxxm)
    nullify(rxym)
    nullify(rxzm)
    nullify(ryym)
    nullify(ryzm)
    nullify(rzzm)
#endif
#endif

  END SUBROUTINE READ_SAVE_FILE_30
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    READHDFMMR 
  !
  ! !DESCRIPTION: read tracer mixing ratio from so-called "old save" file and
  !               uses air mass to convert to mass, in order to initialize the
  !               model.
  !
  !               Used with option istart=4  
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE READHDFMMR( region, file_name, status )
    !
    ! !USES:
    !
    use dims,         only : im, jm, lm, nregions, datadir, region_name
    use io_hdf
    use global_data,  only : mass_dat
    use meteodata  ,  only : m_dat
    use chem_param,   only : fscale, ntracet
    use tm5_distgrid, only : dgrid, Get_DistGrid, scatter
    use partools,     only : isRoot
    !
    ! !INPUT PARAMETERS:
    !
    integer,          intent(in)    :: region
    character(len=*), intent(in)    :: file_name 
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,          intent(out)   :: status
    !
    ! !REVISION HISTORY: 
    !    4 Jul 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition (FIXME: not tested)
    !
    ! !REMARKS:
    ! (1) similar to read_mmix, but all tranported tracers are read at once.
    !     That is: no check on tracers name, and order of tracer matters, 
    !     which depends on the number of processors in TM5 v3.
    ! (2) All tracers must be in the file, and the 4D array should be the
    !     second dataset.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    real, dimension(:,:,:,:), pointer :: rm
    real, dimension(:,:,:),   pointer :: m
    integer                         :: io, sfstart, ifail, sds_id
    integer                         :: n_datasets, n_file_attrs
    integer                         :: istat, attributes, num_type
    integer                         :: sffinfo, sfselect, sfginfo
    integer                         :: sfend, sfrnatt, sfrcatt, sffattr
    integer, dimension(6)           :: idate_save
    integer, dimension(nregions)    :: im_file,jm_file,lm_file
    integer                         :: ntrace_file
    character(len=80)               :: msg_file
    integer                         :: ind,n,i,j,l,offsetn
    
    real,dimension(:,:,:,:), allocatable :: field

    ! start

    m => m_dat(region)%data
    rm => mass_dat(region)%rm
#ifdef slopes
    mass_dat(region)%rxm = 0.
    mass_dat(region)%rym = 0.
    mass_dat(region)%rzm = 0.
#ifdef secmom
    mass_dat(region)%rxxm = 0.
    mass_dat(region)%rxym = 0.
    mass_dat(region)%rxzm = 0.
    mass_dat(region)%ryym = 0.
    mass_dat(region)%ryzm = 0.
    mass_dat(region)%rzzm = 0.
#endif
#endif

    if (isRoot) then 
       allocate(field(im(region),jm(region),lm(region),ntracet))
    else
       allocate(field(1,1,1,1))
    end if


    ! Read on root
    if ( isRoot ) then

       io = sfstart( file_name, DFACC_READ )

       if ( io == -1 ) then
          write(gol,*)'readhdfmmr: Could not open file:'//file_name ; call goErr
          TRACEBACK; status=1; return
       end if
       
       ind = 1  ! select 1st dataset....(0 = m)

       istat = sffinfo(io, n_datasets, n_file_attrs)

       write(gol,*) 'readhdfmmr: Number of datasets and attributes',  &
            n_datasets, n_file_attrs ; call goPr
       
       call io_read4d_32g(io,im(region),jm(region),lm(region),ntracet, &
            field,'rm',ifail,index=ind)
       
       if ( ifail /= 0 .and. ntracet == 1 ) then
          call io_read3d_32g(io,im(region),jm(region),lm(region), &
               field,'rm',ifail,index=ind)
       end if
       
       if ( ifail /= 0 ) then
          write(gol,*) 'readhdfmmr: '//trim(file_name) ; call goErr
          write(gol,*)'readhdfmmr: Failed to read fields'; call goErr
          TRACEBACK; status=1; return
       end if
       write(gol,*) 'readhdfmmr: Read rm, ifail = ', ifail

       istat = sfend(io)
       if (istat /= FAIL) then
          write(gol,*)'readhdfmmr: ... file closed' ; call goErr
          write(gol,*)' ' ; call goPr
       else
          write(gol,*)'readhdfmmr: ERROR in restart from HDF file'; call goErr
          TRACEBACK; status=1; return
       end if

    end if ! root

    
    call SCATTER( dgrid(region), rm, field, mass_dat(region)%halo, status)


    do n=1,ntracet
       rm(:,:,:,n) = rm(:,:,:,n)*m(:,:,:)/fscale(n)
    end do
    
    
    deallocate(field)

    nullify(m)
    nullify(rm)

  END SUBROUTINE READHDFMMR
  !EOC
  
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    READ_MMIX
  !
  ! !DESCRIPTION: read mean mixing ratio of TRANSPORTED tracers from an "mmix hdf
  !               file" and uses air mass to convert to mass. Missing tracers
  !               are set to 0.0
  !               
  !               Used with option istart=5.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE READ_MMIX( region, file_name, status )
    !
    ! !USES:
    !
    use dims,         only : im, jm, lm, nregions, datadir, region_name, parent
    use io_hdf
    use global_data,  only : mass_dat
    use meteodata  ,  only : m_dat
    use chem_param,   only : fscale, ntracet, names
    use tm5_distgrid, only : dgrid, Get_DistGrid, scatter
    use ParTools,     only : isRoot 
    !
    ! !INPUT PARAMETERS:
    !
    integer,          intent(in) :: region 
    character(len=*), intent(in) :: file_name       
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,         intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !    4 Jul 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition (FIXME: not tested)
    !
    ! !REMARKS:
    ! 
    ! (1) This is typically the choice for combining different versions
    !      or extending the number of tracers.
    ! (2) The compounds are searched by name. If not in the mmix file
    !      the field is initialized as zero.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    real,dimension(:,:,:,:),pointer :: rm
#ifdef slopes
    real,dimension(:,:,:,:),pointer :: rxm,rym,rzm
#endif
    real,dimension(:,:,:),pointer   :: m
    integer                         :: io, sfstart, ifail, sds_id
    integer                         :: n_datasets, n_file_attrs
    integer                         :: istat, attributes, num_type
    integer                         :: sffinfo, sfselect, sfginfo
    integer                         :: sfend, sfrnatt, sfrcatt, sffattr
    integer,dimension(6)            :: idate_save
    integer,dimension(nregions)     :: im_file,jm_file,lm_file
    integer                         :: ntrace_file
    character(len=80)               :: msg_file
    integer                         :: ind,n,i,j,l,from_file, my_parent
    real,dimension(:,:,:,:), allocatable :: field
    real,dimension(:,:,:),   allocatable :: field3d

    ! start

    m   =>    m_dat(region)%data
    rm  => mass_dat(region)%rm
#ifdef slopes
    rxm => mass_dat(region)%rxm
    rym => mass_dat(region)%rym
    rzm => mass_dat(region)%rzm
#endif
    
    rm = 0.0
#ifdef slopes
    rxm = 0.0
    rym = 0.0
    rzm = 0.0
#endif

    if (isRoot) then 
       allocate(field(im(region),jm(region),lm(region),ntracet))
    else
       allocate(field(1,1,1,1))
    end if
    

    ! Read on root
    if ( isRoot ) then
       
       from_file = 1

       allocate(field3d(im(region),jm(region),lm(region)))

       io = sfstart(file_name, DFACC_READ)

       if ( io == FAIL ) then 
          my_parent = parent(region)    
          if ( my_parent == 0 ) then
             write(gol,*)'read_mmix : Could not open file and no parent '// &
                  'available :'//file_name ; call goErr
             TRACEBACK; status=1; return
          else

             ! FIXME : this is not going anywhere
             from_file = 0
             write(gol,*) 'readhdf: Trying to initialise '//trim(region_name(region))// &
                  ' from parent....'; call goPr
          endif
       endif
       
       if ( from_file == 1 ) then   ! read from file:

          write(gol,*) 'read_mmix: ',file_name,'... opened for READ access.' ; call goPr

          istat = sffinfo(io, n_datasets, n_file_attrs)

          write(gol,*)  'read_mmix: Number of datasets and attributes', n_datasets, n_file_attrs
          call goPr

          do n = 1, ntracet
             call io_read3d_32(io, im(region), jm(region), lm(region), field3d, names(n), ifail)
             if (ifail == 0) then
                field(1:im(region), 1:jm(region), 1:lm(region), n) = field3d
             else
                write(gol,*) 'Read_mmix:', names(n), '  not found in dataset --> 0.0 ' ; call goPr
                field(1:im(region), 1:jm(region), 1:lm(region), n) = tiny(1.0)
             endif
          enddo

          istat = sfend(io)
          if (istat /= FAIL) then
             write(gol,*)'read_mmix: ... file closed' ; call goPr
             write(gol,*)' '; call goPr
          else
             write(gol,*)'read_mmix: ERROR in restart from HDF file'; call goErr
             TRACEBACK; status=1; return
          end if

          deallocate(field3d) ! no longer needed
       endif  ! from file
    endif ! root   


    call SCATTER( dgrid(region), rm, field, mass_dat(region)%halo, status)


    ! Convert from mmix to kg tracer
    do n = 1, ntracet

       rm(:,:,:,n) = rm(:,:,:,n)*m/fscale(n)

    end do
    
#ifdef slopes
    rxm = 0.0 ; rym = 0.0 ; rzm = 0.0
#endif

    
    ! Done
    deallocate(field)   

    nullify(m, rm)
#ifdef slopes
    nullify(rxm, rym, rzm)
#endif

  END SUBROUTINE READ_MMIX
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    READ_SAVE_FILE
  !
  ! !DESCRIPTION: Read mass of both transpoted AND short lived species from so
  !               -called "save file". No slopes, but regridding to model
  !               resolution is available, since grid information is also read.
  !
  !               Used with istart=31 option.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE READ_SAVE_FILE( region, fname, status, tm4 )
    !
    ! !USES:
    !
    use file_hdf   , only : THdfFile, TSds, Init, Done, ReadAttribute, ReadData
    use grid       , only : TllGridInfo, TLevelInfo, Init, Done, Fill3D
    use dims       , only : im, jm, lm
    use chem_param , only : ntrace, ntracet, fscale, names, ntrace_chem
    use partools   , only : root, ierr, my_real, localComm, isRoot
    use global_data  , only : mass_dat, chem_dat
    use meteodata    , only : global_lli, levi, sp_dat, m_dat
    use tm5_distgrid,  only : dgrid, Get_DistGrid, gather, scatter
    !
    ! !INPUT PARAMETERS:
    !
    integer,           intent(in)           :: region
    character(len=*),  intent(in)           :: fname
    logical,           intent(in), optional :: tm4
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,           intent(out)          ::  status
    !
    ! !REVISION HISTORY: 
    !    4 Oct 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    ! (1) tracers not found in the file are left unchanged (ie to their initialized values)
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    ! 31 (from 60) layer defintition for TM4 save files:
    real, parameter  ::  at_tm4_60_31(0:31) = (/ &
         0.,             7.367743,       467.333588,     1385.912598,    2887.696533,    &
         4941.778320,    6144.314941,    8802.356445,    10209.500977,   11632.758789,   &
         14411.124023,   15706.447266,   16899.468750,   17961.357422,   18864.750000,   &
         19584.330078,   20097.402344,   20384.480469,   20222.205078,   19755.109375,   &
         19027.695313,   18045.183594,   16819.474609,   15379.805664,   13775.325195,   &
         12077.446289,   10376.126953,   8765.053711,    6018.019531,    3960.291504,    &
         2082.273926,    0./)
    real, parameter  ::  bt_tm4_60_31(0:31) = (/ &
         1.00000000,     0.99401945,     0.96764523,     0.93194032,     0.87965691,     &
         0.81125343,     0.77159661,     0.68326861,     0.63554746,     0.58616841,     &
         0.48477158,     0.43396294,     0.38389215,     0.33515489,     0.28832296,     &
         0.24393314,     0.20247590,     0.16438432,     0.09967469,     0.07353383,     &
         0.05169041,     0.03412116,     0.02067788,     0.01114291,     0.00508112,     &
         0.00181516,     0.00046139,     0.00007582,     0.00000000,     0.00000000,     &
         0.00000000,     0.00000000/)

    ! 34 (from 91) layer defintition for TM4 save files:
    real, parameter :: at_tm4_91_34(0:34) = (/ &
         0.000000,     6.575628,    336.772369,   1297.656128,   3010.146973,  &
         5422.802734,  8356.252930, 11543.166992, 14665.645508, 17385.595703, &
         19348.775391, 20319.011719, 20348.916016, 19919.796875, 19184.544922, &
         18191.029297, 16990.623047, 15638.053711, 14192.009766, 12713.897461, &
         11262.484375,  9873.560547,  8564.624023,  7341.469727,  6199.839355,  &
         4663.776367,  3358.425781,  2292.155518,  1463.163940,   857.945801,   &
         450.685791,   204.637451,    76.167656,    21.413612,     0.000000/)

    real, parameter  :: bt_tm4_91_34(0:34) = (/     &
         1.000000,     0.994204,     0.973466,     0.935157,     0.875518,     &
         0.795385,     0.698224,     0.589317,     0.475016,     0.362203,     &
         0.259554,     0.176091,     0.112979,     0.080777,     0.055474,     &
         0.036227,     0.022189,     0.012508,     0.006322,     0.002765,     &
         0.001000,     0.000279,     0.000055,     0.000000,     0.000000,     &
         0.000000,     0.000000,     0.000000,     0.000000,     0.000000,     &
         0.000000,     0.000000,     0.000000,     0.000000,     0.000000/)

    ! 34 (from standard) layer defintition for TM4 save files:
    real, parameter :: at_tm4_std_34(0:34) = (/ &
         0.,              7.367743,        210.393890,      855.361755,     &
         2063.779785,     3850.913330,     6144.314941,     8802.356445,    &
         11632.758789,    14411.124023,    16899.468750,    17961.357422,   &
         18864.750000,    19584.330078,    20097.402344,    20384.480469,   &
         20429.863281,    20222.205078,    19755.109375,    19027.695313,   &
         18045.183594,    16819.474609,    15379.805664,    13775.325195,   &
         12077.446289,    10376.126953,    8765.053711,     7306.631348,    &
         6018.019531,     3960.291504,     1680.640259,     713.218079,     &
         298.495789,      95.636963,       0./)

    real, parameter  :: bt_tm4_std_34(0:34) = (/     &
         1.00000000,      0.99401945,      0.97966272,      0.95182151,     &
         0.90788388,      0.84737492,      0.77159661,      0.68326861,     &
         0.58616841,      0.48477158,      0.38389215,      0.33515489,     &
         0.28832296,      0.24393314,      0.20247590,      0.16438432,     &
         0.13002251,      0.09967469,      0.07353383,      0.05169041,     &
         0.03412116,	  0.02067788,      0.01114291,      0.00508112,     &
         0.00181516,	  0.00046139,      0.00007582,      0.00000000,     &
         0.00000000,      0.00000000,      0.00000000,      0.00000000,     &
         0.00000000,      0.00000000,      0.00000000/)

    ! --- local -------------------------------

    integer                             ::  imr, jmr, lmr

    real, dimension(:,:,:), pointer     ::  sp
    real, dimension(:,:,:), pointer     ::  m

    type(THdfFile)                   ::  hdf
    type(TSds)                       ::  sds_rm, sds_rmc
    integer                          ::  n, ns
    integer                          ::  ntrace_in, ntracet_in
    integer                          ::  n_in, nl
    character(len=8), allocatable    ::  names_in(:)
    integer                          ::  xbeg_in, ybeg_in
    real                             ::  dx_in, dy_in
    integer                          ::  im_in, jm_in, lm_in
    real, allocatable                ::  at_in(:), bt_in(:)
    type(TllGridInfo)                ::  lli_in
    type(TLevelInfo)                 ::  levi_in
    real, allocatable                ::  rm_in(:,:,:,:), sp_gbl(:,:,:)
    real, allocatable                ::  rmt(:,:,:,:),rms(:,:,:,:)
    integer                          ::  l
    logical                          ::  with_tm4

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

    ! global grid sizes
    imr = im(region)
    jmr = jm(region)
    lmr = lm(region)

    ! set pointers
    sp =>   sp_dat(region)%data   ! surface pressure
    m  =>    m_dat(region)%data   ! air mass

    ! special input ?
    with_tm4 = .false.
    if ( present(tm4) ) with_tm4 = tm4


    ! temporary global storage
    if (isRoot) then
       allocate( rmt(imr,jmr,lmr,ntracet    ) )
       allocate( rms(imr,jmr,lmr,ntrace_chem) )
       allocate( sp_gbl(imr,jmr,1) )
    else
       allocate( rmt(1,1,1,1) )
       allocate( rms(1,1,1,1) )
       allocate(sp_gbl(1,1,1))
    endif

    ! init to 0 in case of not found in file
    rmt=tiny(0.)
    rms=tiny(0.)

    
    call gather( dgrid(region), sp, sp_gbl, sp_dat(region)%halo, status) ! get global surf pressure
    IF_NOTOK_RETURN(status=1)

    
    ! read on root (until we implement the stride in reading)
    if ( isRoot ) then

       ! open hdf file
       call Init( hdf, fname, 'read', status )
       IF_NOTOK_RETURN(status=1)

       ! read list with tracer names in file:
       call ReadAttribute( hdf, 'ntrace', ntrace_in, status )
       IF_NOTOK_RETURN(status=1)

       call ReadAttribute( hdf, 'ntracet', ntracet_in, status )
       IF_NOTOK_RETURN(status=1)

       allocate( names_in(ntrace_in) )
       call ReadAttribute( hdf, 'names', names_in, status )
       IF_NOTOK_RETURN(status=1)

       ! read level definition
       call ReadAttribute( hdf, 'lm', lm_in, status )
       IF_NOTOK_RETURN(status=1)

       allocate( at_in(0:lm_in) )
       allocate( bt_in(0:lm_in) )

       if ( with_tm4 ) then
          select case ( lm_in )
          case ( 31 )
             at_in = at_tm4_60_31
             bt_in = bt_tm4_60_31
          case ( 34 )
             at_in = at_tm4_std_34
             bt_in = bt_tm4_std_34
          case default
             write (gol,'("input from TM4 save file implemented for 60-31 and 91-34levels only ...")'); call goErr
             TRACEBACK; status=1; return
          end select
       else
          call ReadAttribute( hdf, 'at', at_in, status )
          IF_NOTOK_RETURN(status=1)
          call ReadAttribute( hdf, 'bt', bt_in, status )
          IF_NOTOK_RETURN(status=1)
       end if

       call Init( levi_in, lm_in, at_in, bt_in, status )
       IF_NOTOK_RETURN(status=1)

       deallocate( at_in )
       deallocate( bt_in )   

       ! read grid definition:
       call ReadAttribute( hdf, 'im', im_in, status )
       IF_NOTOK_RETURN(status=1)
       call ReadAttribute( hdf, 'jm', jm_in, status )
       IF_NOTOK_RETURN(status=1)
       if ( with_tm4 ) then
          xbeg_in = -180.0
          dx_in = 360.0/im_in
          ybeg_in = -90.0
          dy_in = 180.0/jm_in
       else
          call ReadAttribute( hdf, 'xbeg', xbeg_in, status )
          IF_NOTOK_RETURN(status=1)
          call ReadAttribute( hdf, 'dx', dx_in, status )
          IF_NOTOK_RETURN(status=1)
          call ReadAttribute( hdf, 'ybeg', ybeg_in, status )
          IF_NOTOK_RETURN(status=1)
          call ReadAttribute( hdf, 'dy', dy_in, status )
          IF_NOTOK_RETURN(status=1)
       end if
       call Init( lli_in, xbeg_in+0.5*dx_in, dx_in, im_in, &
                          ybeg_in+0.5*dy_in, dy_in, jm_in, status )
       IF_NOTOK_RETURN(status=1)

       ! temporary storage:
       allocate( rm_in(im_in,jm_in,lm_in,1) )


       ! open tracer masses:
       call Init( sds_rm, hdf, 'rm', status )  
       IF_NOTOK_RETURN(status=1)
       if ( .not. with_tm4 ) then
          call Init( sds_rmc, hdf, 'rmc', status )  
          if(status /= 0) write (gol,'("WARNING - No short lived tracer array found in save file...continuing execution")')
          call goPr
       end if

       ! loop over all tracers:
       do n = 1, ntrace

          ! search input record:
          n_in = -1
          do nl = 1, ntrace_in
             if ( names(n) == names_in(nl) ) then
                n_in = nl
                exit
             end if
          end do

          ! not found ? 
          if ( n_in < 0 ) then
             if ( n <= ntracet ) then
                write (gol,'("WARNING - transported tracer not found in save file:")')
             else
                write (gol,'("WARNING - short-lived tracer not found in save file:")')
             end if
             call goPr
             
             write (gol,'("WARNING -   file     : ",a)') trim(fname)              ; call goPr
             write (gol,'("WARNING -   tracer   : ",i3," ",a)') n, trim(names(n)) ; call goPr

             if ( n <= ntracet ) then
                write (gol,'("  initialised tracer ",i2,x,a,"; max mass : ",e10.3)') n, names(n), maxval(rmt(:,:,:,n))
             else
                ns = n-ntracet
                write (gol,'("  initialised tracer ",i2,x,a,"; max mass : ",e10.3)') n, names(n), maxval(rms(:,:,:,ns))
             endif
             call goPr
             cycle
          end if

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

             ! read 3D tracer field from rm or rmc:
             if ( n_in <= ntracet_in ) then
                call ReadData( sds_rm, rm_in, status, start=(/0,0,0,n_in-1/) )  ! kg
                IF_NOTOK_RETURN(status=1)
             else
                call ReadData( sds_rmc, rm_in, status, start=(/0,0,0,n_in-ntracet_in-1/) )  ! kg
                IF_NOTOK_RETURN(status=1)
             end if

             ! expand polar cell for TM4 files ...
             if ( with_tm4 ) then
                do l = 1, lm_in
                   rm_in(:,    1,l,1) = rm_in(1,    1,l,1) / im_in
                   rm_in(:,jm_in,l,1) = rm_in(1,jm_in,l,1) / im_in
                end do
             end if


             ! regrid from saved grid to current model grid
             call Fill3D( global_lli(region), levi, 'n', sp_gbl(1:imr,1:jmr,1), rmt(1:imr,1:jmr,1:lmr,n), &
                          lli_in, levi_in, rm_in(:,:,:,1), 'sum', status )
             IF_NOTOK_RETURN(status=1)

             ! info
             write (gol,'("  initialised tracer ",i2,x,a,"; max mass : ",e10.3)') n, names(n), maxval(rmt(:,:,:,n))
             call goPr
                  
!Need to gather "m" globally to get VMR              
!                  maxval(rmt(1:imr,1:jmr,1:lmr,n)/m_gbl(1:imr,1:jmr,1:lmr)*fscale(n)); call goPr

          else

             ! short lived tracer

             ! read 3D tracer field: 
             if ( with_tm4 ) then
                ! tm4 save files have short lived tracers in rm
                call ReadData( sds_rm, rm_in, status, start=(/0,0,0,n_in-1/) )  ! kg
                IF_NOTOK_RETURN(status=1)
                ! expand polar cell ...
                do l = 1, lm_in
                   rm_in(:,    1,l,1) = rm_in(1,    1,l,1) / im_in
                   rm_in(:,jm_in,l,1) = rm_in(1,jm_in,l,1) / im_in
                end do
             else
                ! read 3D tracer field from rm or rmc:
                if ( n_in <= ntracet_in ) then
                   call ReadData( sds_rm, rm_in, status, start=(/0,0,0,n_in-1/) )  ! kg
                   IF_NOTOK_RETURN(status=1)
                else
                   call ReadData( sds_rmc, rm_in, status, start=(/0,0,0,n_in-ntracet_in-1/) )  ! kg
                   IF_NOTOK_RETURN(status=1)
                end if
             end if

             ! regrid from saved grid to current model grid
             ns = n-ntracet

             call Fill3D( global_lli(region), levi, 'n', sp_gbl(1:imr,1:jmr,1), rms(1:imr,1:jmr,1:lmr,ns), &
                          lli_in,      levi_in, rm_in(:,:,:,1), 'sum', status )
             IF_NOTOK_RETURN(status=1)

             ! info
             write (gol,'("  initialised tracer ",i2,x,a,"; max mass : ",e10.3)') n, names(n), maxval(rms(1:imr,1:jmr,1:lmr,ns))
             call goPr

!Need to gather "m" globally to get VMR              
!               maxval(rms(1:imr,1:jmr,1:lmr,ns)/m(1:imr,1:jmr,1:lmr)*fscale(n)); call goPr

          end if

       end do


       ! cleanup
       deallocate( names_in )
       deallocate( rm_in )

       ! done with grid 
       call Done( lli_in, status ) 
       IF_NOTOK_RETURN(status=1)

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

       ! close data sets
       call Done( sds_rm, status )  
       IF_NOTOK_RETURN(status=1)
       if ( .not. with_tm4 ) then
          call Done( sds_rmc, status )  
       end if

       ! close hdf file (all pe)
       call Done( hdf, status )
       IF_NOTOK_RETURN(status=1)


    end if

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

    call scatter( dgrid(region), chem_dat(region)%rm, rms, chem_dat(region)%halo, status)
    IF_NOTOK_RETURN(status=1)


    ! done
    deallocate( rmt, rms, sp_gbl )
    status = 0

  END SUBROUTINE READ_SAVE_FILE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    WRITE_SAVE_FILE
  !
  ! !DESCRIPTION: save all essential model parameters and fields in a HDF
  !               file. This so-called "save file" is created only if no
  !               restart is written at the end of the run. It can be read
  !               (with READ_SAVE_FILE) to initialize a run, by using
  !               istart=31 option.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE WRITE_SAVE_FILE( msg, file_name, status )
    !
    ! !USES:
    !
    use file_hdf   , only : THdfFile, Init, Done, WriteAttribute, TSds, SetDim, WriteData
    use dims       , only : nregions, region_name
    use dims       , only : im, jm, lm
    use dims       , only : dx, xref, xbeg, xend, ibeg, iend
    use dims       , only : dy, yref, ybeg, yend, jbeg, jend
    use dims       , only : dz, zref, zbeg, zend, lbeg, lend
    use dims       , only : areag
    use dims       , only : czeta, czetak
    use dims       , only : itau, itaui, itaue, itau0, itaut
    use dims       , only : idate, idatei, idatee, idate0, idatet
    use dims       , only : icalendo, iyear0, julday0
    use dims       , only : newyr, newmonth, newday, newsrun
    use dims       , only : tref
    use dims       , only : ndyn, nconv, ndiag, nchem, nsrce
    use dims       , only : nread, ninst, ncheck, ndiff, ndiagp1, ndiagp2, nstep
    use dims       , only : nwrite
    use dims       , only : istart
    use dims       , only : cpu0, cpu1
    use dims       , only : xlabel
    use dims       , only : limits
    use dims       , only : at, bt
    use dims       , only : nsplitsteps, splitorder
    use global_data, only : mass_dat, chem_dat
    use meteodata  , only : m_dat
    use chem_param , only : ntrace, ntrace_chem, ntracet, ra, names, fscale
    use chem_param , only : nstd, istd
    use io_hdf,      only : io_write
    use datetime,    only : tstamp
    use global_data, only : outdir

    use tm5_distgrid,  only : dgrid, Get_DistGrid, GATHER
    use ParTools,      only : isRoot
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*), intent(in)  ::  file_name ! base filename if not found in rc file
    character(len=*), intent(in)  ::  msg       ! message attribute
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          ::  status    ! return code
    !
    ! !REVISION HISTORY: 
    !
    !   11 Nov 2011 - P. Le Sager - adapted for adapted for lon-lat MPI domain decomposition
    !   
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    ! --- local ----------------------------------

    real, dimension(:,:,:,:), pointer       ::  rm
#ifdef slopes
    real, dimension(:,:,:,:), pointer       ::  rxm,rym,rzm
#ifdef secmom
    real, dimension(:,:,:,:), pointer       ::  rxxm,rxym,rxzm,ryym,ryzm,rzzm
#endif
#endif
    real, dimension(:,:,:,:), pointer       ::  rmc
    real, dimension(:,:,:,:), allocatable   ::  glbfield

    type(THdfFile)    ::  hdf
    type(TSds)        ::  sds

    integer           :: region, imr, jmr, lmr, i, n, I1, J1, IH1, halo

    character(len=64)   ::  msg_file
    character(len=256)  ::  fname
    character(len=32 )  ::  key

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

    call InitRc( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    do region = 1,nregions

       if ( isRoot ) then

          ! try if filename for istart.3 was specified in rcfile:

          write (key,'("start.3.",a)') trim(region_name(region))
          call ReadRc( rcF, key, fname, status, default='file_name_empty' )

          if(status == -1) then  ! no specific istart.3 file name found, use default name
             write (fname,'(a,"_",i4.4,3i2.2,"_",a,".hdf")') trim(file_name), idate(1:4), trim(region_name(region))
          else
             write (gol,*) 'Using save file names from rc-file start.3.* values'; call goPr
          endif

          call Init( hdf, trim(fname), 'create', status )
          IF_NOTOK_RETURN(status=1)

          call WriteAttribute( hdf, 'itau'       , itau               , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'nregions'   , nregions           , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'region_name', trim(region_name(region)), status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'im'         , im(region)         , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'jm'         , jm(region)         , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'lm'         , lm(region)         , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'dx'         , dx/xref(region)    , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'dy'         , dy/yref(region)    , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'dz'         , dz/zref(region)    , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'xbeg'       , xbeg(region)       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'xend'       , xend(region)       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'ybeg'       , ybeg(region)       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'yend'       , yend(region)       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'zbeg'       , zbeg(region)       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'zend'       , zend(region)       , status )
          IF_NOTOK_RETURN(status=1)

          if ( region /= 1 ) then
             call WriteAttribute( hdf, 'ibeg', ibeg(region), status )
             IF_NOTOK_RETURN(status=1)
             call WriteAttribute( hdf, 'iend', iend(region), status )
             IF_NOTOK_RETURN(status=1)
             call WriteAttribute( hdf, 'jbeg', jbeg(region), status )
             IF_NOTOK_RETURN(status=1)
             call WriteAttribute( hdf, 'jend', jend(region), status )
             IF_NOTOK_RETURN(status=1)
             call WriteAttribute( hdf, 'lbeg', lbeg(region), status )
             IF_NOTOK_RETURN(status=1)
             call WriteAttribute( hdf, 'lend', lend(region), status )
             IF_NOTOK_RETURN(status=1)
          end if

          call WriteAttribute( hdf, 'xref'       , xref(region), status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'yref'       , yref(region), status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'zref'       , zref(region), status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'tref'       , tref(region), status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'ntracet'    , ntracet     , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'ntrace '    , ntrace      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'nstd'       , nstd        , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'idate'      , idate       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'istart'     , istart      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'ndyn'       , ndyn        , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'nconv'      , nconv       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'ndiag'      , ndiag       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'nchem'      , nchem       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'nsrce'      , nsrce       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'nread'      , nread       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'nwrite'     , nwrite      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'ninst'      , ninst       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'ncheck'     , ncheck      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'ndiff'      , ndiff       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'itaui'      , itaui       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'itaue'      , itaue       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'itaut'      , itaut       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'itau0'      , itau0       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'idatei'     , idatei      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'idatee'     , idatee      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'idatet'     , idatet      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'idate0'     , idate0      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'icalendo'   , icalendo    , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'iyear0'     , iyear0      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'julday0'    , julday0     , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'ndiagp1'    , ndiagp1     , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'ndiagp2'    , ndiagp2     , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'nstep'      , nstep       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'cpu0'       , cpu0        , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'cpu1'       , cpu1        , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'ra'         , ra          , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'fscale'     , fscale      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'names'      , names       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'areag'      , areag       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'czeta'      , czeta       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'czetak'     , czetak      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'xlabel'     , xlabel      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'istd'       , istd        , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'newyr'      , newyr       , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'newmonth'   , newmonth    , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'newday'     , newday      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'newsrun'    , newsrun     , status )
          IF_NOTOK_RETURN(status=1)
!          call WriteAttribute( hdf, 'cdebug'     , cdebug      , status )
!          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'limits'     , limits      , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'at'         , at          , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'bt'         , bt          , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'nsplitsteps', nsplitsteps , status )
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( hdf, 'splitorder' , splitorder  , status )
          IF_NOTOK_RETURN(status=1)

          msg_file = msg
          call WriteAttribute( hdf, 'msg', msg_file, status )
          IF_NOTOK_RETURN(status=1)

       end if  ! root all PEs from here


       ! Short cuts
       !----------------------
       rm =>   mass_dat(region)%rm
#ifdef slopes
       rxm =>  mass_dat(region)%rxm
       rym =>  mass_dat(region)%rym
       rzm =>  mass_dat(region)%rzm
#ifdef secmom
       rxxm => mass_dat(region)%rxxm
       rxym => mass_dat(region)%rxym
       rxzm => mass_dat(region)%rxzm
       ryym => mass_dat(region)%ryym
       ryzm => mass_dat(region)%ryzm
       rzzm => mass_dat(region)%rzzm
#endif
#endif

       ! Local & global bounds
       !-----------------------
       call Get_DistGrid( dgrid(region), I_STRT=i1, J_STRT=j1, I_STRT_HALO=ih1 )
       halo = i1-ih1
       imr = im(region) ; jmr = jm(region) ; lmr = lm(region)


       ! Allocate global array
       !----------------------
       if (isRoot) then 
          allocate( glbfield( imr, jmr, lmr, ntracet) )
          glbfield = 0.0
       else
          allocate( glbfield(1,1,1,1) )
       end if


       ! M (Air mass)
       !----------------------
       CALL GATHER( dgrid(region), m_dat(region)%data, glbfield(:,:,:,1), m_dat(region)%halo, status )
       IF_ERROR_RETURN(status=1)

       if ( isRoot ) then 

          call Init( sds, hdf, 'm', (/imr,jmr,lmr/), 'real(4)', status )
          call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status)
          call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status)
          call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status)
          call WriteData(sds, glbfield(:,:,:,1), status)
          call Done(sds,status)
          IF_ERROR_RETURN(status=1)

       end if


       ! RM (transported tracer masses)
       !-----------------------------
       CALL GATHER( dgrid(region), rm, glbfield, halo, status )
       IF_ERROR_RETURN(status=1)

       if ( isRoot ) then 

          call Init( sds, hdf, 'rm', (/imr,jmr,lmr,ntracet/), 'real(4)', status )
          call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status)
          call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status)
          call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status)
          call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status)
          call WriteData(sds,glbfield,status)
          call Done(sds,status)
          IF_ERROR_RETURN(status=1)

       end if


#ifdef slopes

       ! RXM
       !-----------------------------
       CALL GATHER( dgrid(region), rxm, glbfield, halo, status )
       IF_ERROR_RETURN(status=1)

       if ( isRoot ) then

          call Init( sds, hdf, 'rxm', (/imr,jmr,lmr,ntracet/), 'real(4)', status )
          call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status)
          call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status)
          call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status)
          call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status)
          call WriteData(sds,glbfield,status)
          call Done(sds,status)
          IF_ERROR_RETURN(status=1)

       end if

       ! RYM
       !-----------------------------
       CALL GATHER( dgrid(region), rym, glbfield, halo, status )
       IF_ERROR_RETURN(status=1)

       if ( isRoot ) then

          call Init( sds, hdf, 'rym', (/imr,jmr,lmr,ntracet/), 'real(4)', status )
          call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status)
          call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status)
          call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status)
          call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status)
          call WriteData(sds,glbfield,status)
          call Done(sds,status)
          IF_ERROR_RETURN(status=1)

       end if

       ! RZM
       !-----------------------------
       CALL GATHER( dgrid(region), rzm, glbfield, halo, status )
       IF_ERROR_RETURN(status=1)

       if ( isRoot ) then

          call Init( sds, hdf, 'rzm', (/imr,jmr,lmr,ntracet/), 'real(4)', status )
          call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status)
          call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status)
          call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status)
          call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status)
          call WriteData(sds,glbfield,status)
          call Done(sds,status)
          IF_ERROR_RETURN(status=1)

       end if

#ifdef secmom
       ! second moments

       ! RXXM
       !-----------------------------
       CALL GATHER( dgrid(region), rxxm, glbfield, halo, status )
       IF_ERROR_RETURN(status=1)

       if ( isRoot ) then
          call Init( sds, hdf, 'rxxm', (/imr,jmr,lmr,ntracet/), 'real(4)', status )
          call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status)
          call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status)
          call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status)
          call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status)
          call WriteData(sds,glbfield,status)
          call Done(sds,status)
          IF_ERROR_RETURN(status=1)
       end if

       ! RXYM
       !-----------------------------
       CALL GATHER( dgrid(region), rxym, glbfield, halo, status )
       IF_ERROR_RETURN(status=1)

       if ( isRoot ) then
          call Init( sds, hdf, 'rxym', (/imr,jmr,lmr,ntracet/), 'real(4)', status )
          call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status)
          call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status)
          call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status)
          call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status)
          call WriteData(sds,glbfield,status)
          call Done(sds,status)
          IF_ERROR_RETURN(status=1)
       end if

       ! RXZM
       !-----------------------------
       CALL GATHER( dgrid(region), rxzm, glbfield, halo, status )
       IF_ERROR_RETURN(status=1)

       if ( isRoot ) then
          call Init( sds, hdf, 'rxzm', (/imr,jmr,lmr,ntracet/), 'real(4)', status )
          call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status)
          call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status)
          call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status)
          call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status)
          call WriteData(sds,glbfield,status)
          call Done(sds,status)
          IF_ERROR_RETURN(status=1)
       end if

       ! RYYM
       !-----------------------------
       CALL GATHER( dgrid(region), ryym, glbfield, halo, status )
       IF_ERROR_RETURN(status=1)

       if ( isRoot) then
          call Init( sds, hdf, 'ryym', (/imr,jmr,lmr,ntracet/), 'real(4)', status )
          call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status)
          call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status)
          call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status)
          call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status)
          call WriteData(sds,glbfield,status)
          call Done(sds,status)
          IF_ERROR_RETURN(status=1)
       end if

       ! RYZM
       !-----------------------------
       CALL GATHER( dgrid(region), ryzm, glbfield, halo, status )
       IF_ERROR_RETURN(status=1)

       if ( isRoot) then
          call Init( sds, hdf, 'ryzm', (/imr,jmr,lmr,ntracet/), 'real(4)', status )
          call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status)
          call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status)
          call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status)
          call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status)
          call WriteData(sds,glbfield,status)
          call Done(sds,status)
          IF_ERROR_RETURN(status=1)
       endif

       ! RZZM
       !-----------------------------
       CALL GATHER( dgrid(region), rzzm, glbfield, halo, status )
       IF_ERROR_RETURN(status=1)

       if ( isRoot) then
          call Init( sds, hdf, 'rzzm', (/imr,jmr,lmr,ntracet/), 'real(4)', status )
          call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status)
          call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status)
          call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status)
          call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status)
          call WriteData(sds,glbfield,status)
          call Done(sds,status)
          IF_ERROR_RETURN(status=1)
       endif

#endif
#endif

       ! clear
       deallocate(glbfield)
       nullify(rm)
#ifdef slopes
       nullify(rxm)
       nullify(rym)
       nullify(rzm)
#ifdef secmom
       nullify(rxxm)
       nullify(rxym)
       nullify(rxzm)
       nullify(ryym)
       nullify(ryzm)
       nullify(rzzm)
#endif
#endif

       ! SHORT LIVED TRACERS if any
       ! ----------------------------
       if ( ntrace_chem > 0 ) then

          ! allocate global array
          if (isRoot) then 
             allocate(glbfield(imr,jmr,lmr,ntrace_chem))
          else
             allocate(glbfield(1,1,1,1))
          end if

          ! tracer mass
          rmc => chem_dat(region)%rm

          CALL GATHER( dgrid(region), rmc, glbfield, 0, status )
          IF_NOTOK_RETURN(status=1)

          if ( isRoot ) then

             ! write tracers
             call io_write( hdf%id, imr,'LON'//region_name(region),&
                            jmr,'LAT'//region_name(region), &
                            lmr,'HYBRID',ntrace_chem,'NTRACE_CHEM',glbfield,'rmc')

          end if

          ! clear
          nullify(rmc)      
          deallocate(glbfield)

       end if  ! ntrace_chem > 0


       ! Close
       !-------
       if ( isRoot ) then
          call Done( hdf, status )
          IF_NOTOK_RETURN(status=1)
       end if

    end do  !regions...


    ! -------
    ! Done
    ! -------
    call DoneRc( rcF, status )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0

  END SUBROUTINE WRITE_SAVE_FILE
  !EOC

END MODULE IO_SAVE