!#################################################################
!
! contains routines to read emissions and 
! to read and write the main model state from/to file
!
!### 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
!
#include "tm5.inc"
!
!#################################################################

module io_save

  use GO, only : gol, goErr, goPr
  use GO, only : TrcFile, InitRc => Init, DoneRc => Done, ReadRc
  use global_data,   only : rcfile
  
  implicit none
  
  ! --- in/out --------------------------

  private

  public :: readhdf, readhdfmmr, savehdf, readhdfmmix
  public :: readtm3hdf, readtm3hdf_4d


  ! --- const ----------------------------------

  character(len=*), parameter  ::  mname = 'io_save'
  type(TrcFile)       ::  rcF


contains


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

  !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
          print *, 'readtm3hdf: io_read2d_32g failed:'
          print *, 'readtm3hdf:   filename   : ', trim(filename)
          print *, 'readtm3hdf:   rank       : ', rank
          print *, 'readtm3hdf:   im_emis    : ', im_emis
          print *, 'readtm3hdf:   jm_emis    : ', jm_emis
          print *, 'readtm3hdf:   lm_emis    : ', lm_emis
          print *, 'readtm3hdf:   month      : ', month
          print *, 'readtm3hdf:   field_name : ', trim(field_name)
          if (present(idt)) print *, 'readtm3hdf:   idt        : ', idt
          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
          print *, 'readtm3hdf: io_read2d_32g failed:'
          print *, 'readtm3hdf:   filename   : ', trim(filename)
          print *, 'readtm3hdf:   rank       : ', rank
          print *, 'readtm3hdf:   im_emis    : ', im_emis
          print *, 'readtm3hdf:   jm_emis    : ', jm_emis
          print *, 'readtm3hdf:   lm_emis    : ', lm_emis
          print *, 'readtm3hdf:   month      : ', month
          print *, 'readtm3hdf:   field_name : ', trim(field_name)
          if (present(idt)) print *, 'readtm3hdf:   idt        : ', idt
          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
          print *, 'readtm3hdf: io_read2d_32g failed:'
          print *, 'readtm3hdf:   filename   : ', trim(filename)
          print *, 'readtm3hdf:   rank       : ', rank
          print *, 'readtm3hdf:   im_emis    : ', im_emis
          print *, 'readtm3hdf:   jm_emis    : ', jm_emis
          print *, 'readtm3hdf:   lm_emis    : ', lm_emis
          print *, 'readtm3hdf:   month      : ', month
          print *, 'readtm3hdf:   field_name : ', trim(field_name)
          if (present(idt)) print *, 'readtm3hdf:   idt        : ', idt
          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!
           print *, '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
             print *, 'readtm3hdf: io_read2d_32g failed:'
             print *, 'readtm3hdf:   filename   : ', trim(filename)
             print *, 'readtm3hdf:   rank       : ', rank
             print *, 'readtm3hdf:   im_emis    : ', im_emis
             print *, 'readtm3hdf:   jm_emis    : ', jm_emis
             print *, 'readtm3hdf:   lm_emis    : ', lm_emis
             print *, 'readtm3hdf:   month      : ', month
             print *, 'readtm3hdf:   field_name : ', trim(field_name)
             if (present(idt)) print *, 'readtm3hdf:   idt        : ', idt
             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

  ! ***

  subroutine readtm3hdf_4d( filename, rank, im_emis, jm_emis, lm_emis, km_emis, &
                         month, field, field_name, idt )
                         
    use file_hdf, only : THdfFile, Init, Done
    use io_hdf  , only : FAIL
    use io_hdf  , only : io_read4d_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,km_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(:,:,:)
    real, allocatable        ::  field4d(:,:,:,:)

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

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

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

    select case(rank) 

    case (4)   ! 4D field

       allocate(field4d(im_emis,jm_emis,lm_emis,km_emis))
       if(present(idt)) call io_read4d_32g(hdf%id,im_emis,jm_emis,lm_emis,km_emis,field4d,field_name, status,idate=idt)
       if(.not.present(idt)) call io_read4d_32g(hdf%id,im_emis,jm_emis,lm_emis,km_emis,field4d,field_name, status,index=idx)
       if ( status /= 0 )  then
          print *, 'readtm3hdf_4d: io_read4d_32g failed:'
          print *, 'readtm3hdf_4d:   filename   : ', trim(filename)
          print *, 'readtm3hdf_4d:   rank       : ', rank
          print *, 'readtm3hdf_4d:   im_emis    : ', im_emis
          print *, 'readtm3hdf_4d:   jm_emis    : ', jm_emis
          print *, 'readtm3hdf_4d:   lm_emis    : ', lm_emis
          print *, 'readtm3hdf_4d:   month      : ', month
          print *, 'readtm3hdf_4d:   field_name : ', trim(field_name)
          if (present(idt)) print *, 'readtm3hdf_4d:   idt        : ', idt
          call escape_tm('readtm3hdf_4d: Error reading HDF file 3D')
       endif
       field(:,:,:)=field4d(:,:,:,month)
       deallocate(field4d)

    case default
       write (gol,*) 'readtm3hdf_4d: Please check the rank of the required field,', &
               ' must be 4'; call goPr
    end select

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

  end subroutine readtm3hdf_4d

  ! ***
  

  !
  ! data are expected in kg   
  ! (typical save.hdf files from previous run...)
  ! only long-lived tracer are read, since the fields 
  ! for chemistry are not yet allocated
  !

  subroutine ReadHDF( region, file_name )

    use io_hdf
    use dims,        only : region_name, im, jm, lm, nregions, &
                            datadir, parent
    use global_data, only : mass_dat
    use meteodata  , only : m_dat
    use chem_param,  only : ntrace, ntracet
    use toolbox,     only : escape_tm
    use ParTools   , only : myid, root_t
    use ParTools   , only : ntracet_ar, ntracetloc
#ifdef MPI
    use mpi_comm,    only : scatter_after_read_t
    use mpi_const,   only : com_trac, ierr, my_real, mpi_integer
#endif

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

    integer, intent(in)             :: region
    character(len=*), intent(in)    :: file_name 
    
    ! --- local -------------------------------

    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, 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, i,j,l,n
    character(len=80)               :: msg_file
    integer                         :: ind, my_parent,imr,jmr,lmr
    real,dimension(:,:,:,:), allocatable :: field, field2
    real,dimension(:,:,:), allocatable   :: msave
    integer                         :: from_file  ! signals reading from file

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

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

#ifdef slopes
    rxm => mass_dat(region)%rxm_t
    rym => mass_dat(region)%rym_t
    rzm => mass_dat(region)%rzm_t
#ifdef secmom
    rxxm => mass_dat(region)%rxxm_t
    rxym => mass_dat(region)%rxym_t
    rxzm => mass_dat(region)%rxzm_t
    ryym => mass_dat(region)%ryym_t
    ryzm => mass_dat(region)%ryzm_t
    rzzm => mass_dat(region)%rzzm_t
#endif
#endif


    imr = im(region) ; jmr = jm(region) ; lmr = lm(region)

    allocate(field(-1:imr+2,-1:jmr+2,lmr,ntracet))
    field = 0.0

    allocate(msave(imr,jmr,lmr))

    if ( myid == root_t ) then

       from_file = 1
       allocate(field2(imr,jmr,lmr,ntracet))   
       ind = len_trim(datadir)
       !io = sfstart(datadir(1:ind)//file_name//trim(region_name(region)), DFACC_READ)
       io = sfstart( file_name, DFACC_READ )
       if ( io == FAIL ) then 
          my_parent = parent(region)    
          if ( my_parent == 0 ) then
             !call escape_tm('readhdf: Could not open file and no parent '// &
             !     'available :'//datadir(1:ind)//file_name//trim(region_name(region)))
             call escape_tm('readhdf: Could not open file and no parent '// &
                  'available :'//trim(file_name) )
          else
             from_file = 0
             write (*,'("readhdf: Trying to initialise ",a," from parent....")') trim(region_name(region))
          end if
       end if
       if ( from_file == 1 ) then
          write (*,'("readhdf:")')
          write (*,'("readhdf: opened for READ access: ",a)') trim(file_name)

          call io_read3d_32g(io,imr,jmr,lmr,msave,'m',ifail)
          if ( ifail /= 0 ) then
             call escape_tm('readhdf: Failed to read m from saveold.hdf in readhdf')
          end if

          call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rm',ifail)
          if( ifail /= 0 .and. ntracet == 1 ) then
             call io_read3d_32g(io,imr,jmr,lmr,field2,'rm',ifail)
          end if

          if ( ifail /= 0 ) &
               call escape_tm('readhdf: Failed to read in saveold.hdf in readhdf')
          write (*,'("readhdf: Read rm, ifail = ",i6)') ifail

          field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers!
       end if  !from_file

    end if  !root_t

#ifdef MPI
    ! broadcast from_file !
    call mpi_bcast(from_file, 1, mpi_integer , root_t, com_trac, ierr)
#endif

    if ( from_file == 1 ) then 
#ifdef MPI   
       call mpi_bcast(msave, imr*jmr*lmr, my_real, root_t, com_trac, ierr)
       call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rm,root_t)
#else
       rm = field
#endif
    else
       call init_child(region)   ! get info from parent....
    end if

    field = 0.0

    if ( from_file == 1 ) then   
    
#ifdef slopes

          if ( myid == root_t ) then

             call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rxm',ifail)
             if ( ifail /= 0 .and. ntracet == 1 ) then
                call io_read3d_32g(io,imr,jmr,lmr,field2,'rxm',ifail)
             end if
             if ( ifail /= 0 ) then
                call escape_tm( &
                  'readhdf: Failed to read in rxm in saveold.hdf in readhdf')
             end if
             write (*,'("readhdf: Read rxm, ifail = ",i6)') ifail

             field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers!

          end if  !root_t

#ifdef MPI   
          call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rxm,root_t)
#else
          rxm = field
#endif
          field = 0.0

          if(myid == root_t) then

             call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rym',ifail)
             if ( ifail /= 0 .and. ntracet == 1 ) then
                call io_read3d_32g(io,imr,jmr,lmr,field2,'rym',ifail)
             end if
             if ( ifail /= 0 ) then
                call escape_tm( &
                     'readhdf: Failed to read in rym in saveold.hdf in readhdf')
             end if
             write (*,'("readhdf: Read rym, ifail = ",i6)') ifail
             field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers!

          end if  !root_t

#ifdef MPI   
          call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rym,root_t)
#else
          rym = field
#endif
          field = 0.0

          if ( myid == root_t ) then

             call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rzm',ifail)
             if ( ifail /= 0 .and. ntracet == 1 ) then
                call io_read3d_32g(io,imr,jmr,lmr,field2,'rzm',ifail)
             end if
             if ( ifail /= 0 ) then
                call escape_tm( &
                     'readhdf: Failed to read in rzm in saveold.hdf in readhdf')
             end if
             write (*,'("readhdf: Read rzm, ifail = ",i6)') ifail

             field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers!
             deallocate(field2)

          end if  !root_t

#ifdef MPI   
          call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rzm,root_t)
#else
          rzm = field
#endif

#endif


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

#ifdef secmom
         allocate(field2(imr,jmr,lmr,ntracet))   

       if ( myid == root_t) then

         call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rxxm',ifail)
         if ( ifail /= 0 .and. ntracet == 1 ) then
           call io_read3d_32g(io,imr,jmr,lmr,field2,'rxxm',ifail)
         endif
         if ( ifail /= 0 ) then
           call escape_tm('readhdf: Failed to read in rxxm in saveold.hdf in readhdf')
         end if
         write (*,'("readhdf: Read rxxm, ifail = ",i6)') ifail

        field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers!

       endif  !root_t

#ifdef MPI   
   call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rxxm,root_t)
#else
   rxxm = field
#endif
   field = 0.0

   if(myid == root_t) then

      call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rxym',ifail)
             if ( ifail /= 0 .and. ntracet == 1 ) then
         call io_read3d_32g(io,imr,jmr,lmr,field2,'rxym',ifail)
      endif
             if ( ifail /= 0 ) then
                call escape_tm( &
                     'readhdf: Failed to read in rxym in saveold.hdf in readhdf')
             end if
      write (*,'("readhdf: Read rxym, ifail = ",i6)') ifail
      field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers!

   endif  !root_t

#ifdef MPI   
   call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rxym,root_t)
#else
   rxym = field
#endif
   field = 0.0
 
   if(myid == root_t) then

      call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rxzm',ifail)
             if ( ifail /= 0 .and. ntracet == 1 ) then
         call io_read3d_32g(io,imr,jmr,lmr,field2,'rxzm',ifail)
      endif
             if ( ifail /= 0 ) then
                call escape_tm( &
                     'readhdf: Failed to read in rxzm in saveold.hdf in readhdf')
             end if
             write (*,'("readhdf: Read rxzm, ifail = ",i6)') ifail

      field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers!

   endif  !root_t

#ifdef MPI   
   call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rxzm,root_t)
#else
   rxzm = field
#endif

   if(myid == root_t) then

      call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'ryym',ifail)
             if ( ifail /= 0 .and. ntracet == 1 ) then
         call io_read3d_32g(io,imr,jmr,lmr,field2,'ryym',ifail)
      endif
             if ( ifail /= 0 ) then
                call escape_tm( &
                     'readhdf: Failed to read in ryym in saveold.hdf in readhdf')
             end if
      write (*,'("readhdf: Read ryym, ifail = ",i6)') ifail

      field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers!

   endif  !root_t

#ifdef MPI   
   call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,ryym,root_t)
#else
   ryym = field
#endif
   field = 0.0

   if(myid == root_t) then

      call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'ryzm',ifail)
             if ( ifail /= 0 .and. ntracet == 1 ) then
         call io_read3d_32g(io,imr,jmr,lmr,field2,'ryzm',ifail)
      endif
             if ( ifail /= 0 ) then
                call escape_tm( &
                     'readhdf: Failed to read in ryzm in saveold.hdf in readhdf')
             end if
      write (*,'("readhdf: Read ryzm, ifail = ",i6)') ifail
      field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers!

   endif  !root_t

#ifdef MPI   
   call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,ryzm,root_t)
#else
   ryzm = field
#endif
   field = 0.0
 
   if(myid == root_t) then

      call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rzzm',ifail)
             if ( ifail /= 0 .and. ntracet == 1 ) then
         call io_read3d_32g(io,imr,jmr,lmr,field2,'rzzm',ifail)
      endif
             if ( ifail /= 0 ) then
                call escape_tm( &
                     'readhdf: Failed to read in rzzm in saveold.hdf in readhdf')
             end if
          write (*,'("readhdf: Read rzzm, ifail = ",i6)') ifail

      field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers!
      deallocate(field2)

   endif  !root_t

#ifdef MPI   
   call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rzzm,root_t)
#else
   rzzm = field
#endif
 
#endif

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

       if ( myid == root_t ) then

          istat = sfend(io)
          if ( istat /= FAIL ) then
             write (*,'("readhdf: ... file closed")')
             write (*,'("readhdf: ")')
          else
             call escape_tm('readhdf: ERROR in restart from HDF file')
          end if

       end if

       ! scale in input by a possible factor.
       ! rm/msave ---> mixing ratio on file. mr*m --> something back in kilo's

       if ( myid == root_t ) then

          print *, 'readhdf: Maximum ratio of the air masses:', &
               maxval(abs(msave/m(1:imr,1:jmr,1:lmr)))
          print *, 'readhdf: Minimum ratio of the air masses:', &
               minval(abs(msave/m(1:imr,1:jmr,1:lmr)))

       end if

       do n=1,ntracetloc

          do l = 1,lmr
             do j = 1,jmr
                do i = 1,imr
                   rm(i,j,l,n)  = rm (i,j,l,n)*m(i,j,l)/msave(i,j,l)
#ifdef slopes
                   rxm(i,j,l,n) = rxm(i,j,l,n)*m(i,j,l)/msave(i,j,l)
                   rym(i,j,l,n) = rym(i,j,l,n)*m(i,j,l)/msave(i,j,l)
                   rzm(i,j,l,n) = rzm(i,j,l,n)*m(i,j,l)/msave(i,j,l)
#ifdef secmom
                   rxxm(i,j,l,n) = rxxm(i,j,l,n)*m(i,j,l)/msave(i,j,l)
                   rxym(i,j,l,n) = rxym(i,j,l,n)*m(i,j,l)/msave(i,j,l)
                   rxzm(i,j,l,n) = rxzm(i,j,l,n)*m(i,j,l)/msave(i,j,l)
                   ryym(i,j,l,n) = ryym(i,j,l,n)*m(i,j,l)/msave(i,j,l)
                   ryzm(i,j,l,n) = ryzm(i,j,l,n)*m(i,j,l)/msave(i,j,l)
                   rzzm(i,j,l,n) = rzzm(i,j,l,n)*m(i,j,l)/msave(i,j,l)
#endif
#endif
                end do
             end do
          end do
       end do
       
    end if   ! from_file == 1

    deallocate(msave)
    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 ReadHDF


  ! ***
  

  subroutine init_child(region)
    !
    ! Parent is used to initilise rm, rxm,rym,rzm of region. 
    ! Assumed parallel over tracers...
    !
    use dims,        only : lm, ibeg, iend, jbeg, jend, xref, yref, parent
    use tracer_data, only : mass_dat
    use meteodata  , only : m_dat
    use chem_param,  only : ntrace, ntracet
    use toolbox,     only : escape_tm
    use ParTools,    only : myid, ntracetloc, root_t

    implicit none

    ! in/out
    integer, intent(in) :: region

    ! local
    real,dimension(:,:,:,:),pointer :: rm
#ifdef slopes
    real,dimension(:,:,:,:),pointer :: rxm,rym,rzm
#endif
    real,dimension(:,:,:),pointer   :: m
    real,dimension(:,:,:,:),pointer :: rmp   !parent
    real,dimension(:,:,:),pointer   :: mp

    real           :: mm_parent
    integer        :: my_parent
    integer        :: n,l,jp,ip,ic,jc, xref_, yref_, i,j

    ! start

    m => m_dat(region)%data
    rm => mass_dat(region)%rm_t
#ifdef slopes
    rxm => mass_dat(region)%rxm_t
    rym => mass_dat(region)%rym_t
    rzm => mass_dat(region)%rzm_t
#endif

    my_parent = parent(region)    

    if ( my_parent == 0 ) then
       call escape_tm( 'init_child: Region 1 cannot be initialised')
    end if

    mp => m_dat(my_parent)%data
    rmp => mass_dat(my_parent)%rm_t

    xref_ = xref(region)/xref(my_parent)
    yref_ = yref(region)/yref(my_parent)
    do n=1,ntracetloc
       do l=1,lm(region)
          do jp=jbeg(region), jend(region)
             jc = (jp-jbeg(region))*yref_
             do ip=ibeg(region), iend(region)
                mm_parent = rmp(ip,jp,l,n)/mp(ip,jp,l)
                ic = (ip-ibeg(region))*xref_
                do j=1,yref_
                   do i=1,xref_
                      rm (ic+i,jc+j,l,n) = mm_parent*m(ic+i,jc+j,l)  
#ifdef slopes
                      rxm(ic+i,jc+j,l,n) = 0.0
                      rym(ic+i,jc+j,l,n) = 0.0
                      rzm(ic+i,jc+j,l,n) = 0.0
#endif
                   end do
                end do
             end do
          end do
       end do
    end do
    if ( myid == root_t ) print *, 'readhdf: Initialisation successful '

    nullify(m)
    nullify(rm)
#ifdef slopes
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
#endif
    nullify(mp)
    nullify(rmp)
    
  end subroutine init_child

  subroutine init_child_edge(region)
    !
    ! ONLY the edges!
    ! Parent is used to initilise rm, rxm,rym,rzm of region. 
    ! Assumed parallel over tracers...
    !
    use dims,        only : lm, ibeg, iend, jbeg, jend, xref, yref, parent, im, jm
    use global_data, only : region_dat
    use tracer_data, only : mass_dat
    use meteodata  , only : m_dat
    use chem_param,  only : ntrace, ntracet
    use toolbox,     only : escape_tm
    use ParTools,    only : myid, ntracetloc, root_t

    implicit none

    ! in/out
    integer, intent(in) :: region

    ! local
    real,dimension(:,:,:,:),pointer :: rm
#ifdef slopes
    real,dimension(:,:,:,:),pointer :: rxm,rym,rzm
#endif
    real,dimension(:,:,:),pointer   :: m
    real,dimension(:,:,:,:),pointer :: rmp   !parent
    real,dimension(:,:,:),pointer   :: mp
    integer,dimension(:,:),pointer  :: edge

    real           :: mm_parent
    integer        :: my_parent
    integer        :: n,l,jp,ip,ic,jc, xref_, yref_, i,j

    ! start

    my_parent = parent(region)    
    if ( my_parent == 0 ) return

    m => m_dat(region)%data
    rm => mass_dat(region)%rm_t
#ifdef slopes
    rxm => mass_dat(region)%rxm_t
    rym => mass_dat(region)%rym_t
    rzm => mass_dat(region)%rzm_t
#endif
    edge => region_dat(region)%edge


    mp => m_dat(my_parent)%data
    rmp => mass_dat(my_parent)%rm_t

    xref_ = xref(region)/xref(my_parent)
    yref_ = yref(region)/yref(my_parent)
    do n=1,ntracetloc
       do l=1,lm(region)
          do j=1,jm(region)
             jp = jbeg(region) + (j-1)/yref_
             do i=1,im(region)
                if (edge(i,j) /= -1 ) cycle
                ip = ibeg(region) + (i-1)/xref_
                mm_parent = rmp(ip,jp,l,n)/mp(ip,jp,l)
                rm (i,j,l,n) = mm_parent*m(i,j,l)  
#ifdef slopes
                rxm(i,j,l,n) = 0.0
                rym(i,j,l,n) = 0.0
                rzm(i,j,l,n) = 0.0
#endif
             end do
          end do
       end do
    end do
    if ( myid == root_t ) print *, 'readhdf: Edges Initialisation successful '

    nullify(m)
    nullify(rm)
#ifdef slopes
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
#endif
    nullify(mp)
    nullify(rmp)
    
  end subroutine init_child_edge





  subroutine readhdfmmr(region,file_name)
    !
    ! data are expected in mixing ratio...
    !
    use dims,        only : im, jm, lm, nregions, datadir, region_name
    use io_hdf
    use tracer_data, only : mass_dat
    use meteodata  , only : m_dat
    use chem_param,  only : fscale, ntracet
    use toolbox,     only : escape_tm
#ifdef MPI
    use mpi_comm,    only : scatter_after_read_t
    use mpi_const,   only : root, myid, ntracet_ar, ntracetloc
#endif

    implicit none

    ! in/out
    integer, intent(in)             :: region
    character(len=*), intent(in)    :: file_name 

    ! 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   :: 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,field2

    ! start

    m => m_dat(region)%data
    rm => mass_dat(region)%rm_t
#ifdef slopes
    rxm => mass_dat(region)%rxm_t
    rym => mass_dat(region)%rym_t
    rzm => mass_dat(region)%rzm_t
#ifdef secmom
    rxxm => mass_dat(region)%rxxm_t
    rxym => mass_dat(region)%rxym_t
    rxzm => mass_dat(region)%rxzm_t
    ryym => mass_dat(region)%ryym_t
    ryzm => mass_dat(region)%ryzm_t
    rzzm => mass_dat(region)%rzzm_t
#endif
#endif
    allocate(field(-1:im(region)+2,-1:jm(region)+2,1:lm(region),1:ntracet))
    field=0.0
#ifdef MPI
    if ( myid == root ) then 
#endif
       ! to hold all ntrace values of rm
       allocate(field2(im(region),jm(region),lm(region),ntracet))
       io = sfstart( file_name, DFACC_READ )

       if ( io == -1 ) call escape_tm('readhdfmmr: Could not open file:'//file_name)
       ind = 1  ! select 1st dataset....(0 = m)

       istat = sffinfo(io, n_datasets, n_file_attrs)
       print*, 'readhdfmmr: Number of datasets and attributes',  &
            n_datasets, n_file_attrs
       call io_read4d_32g(io,im(region),jm(region),lm(region),ntracet, &
            field2,'rm',ifail,index=ind)
       if ( ifail /= 0 .and. ntracet == 1 ) then
          call io_read3d_32g(io,im(region),jm(region),lm(region), &
               field2,'rm',ifail,index=ind)
       end if
       if ( ifail /= 0 ) then
          print *, 'readhdfmmr: '//trim(file_name)
          call escape_tm('readhdfmmr: Failed to read fields')
       end if
       print*, 'readhdfmmr: Read rm, ifail = ', ifail

       istat = sfend(io)
       if (istat /= FAIL) then
          print*,'readhdfmmr: ... file closed'
          print*,' '
       else
          call escape_tm('readhdfmmr: ERROR in restart from HDF file')
       end if

       ! only transported tracers!
       field(1:im(region),1:jm(region),1:lm(region),:)=field2
       deallocate(field2) ! no longer needed
#ifdef MPI
    end if ! root
    call scatter_after_read_t(field,im(region),jm(region),lm(region),&
         2,2,0,rm,root)
#else
    rm = field
#endif
#ifdef MPI
    offsetn = sum(ntracet_ar(0:myid-1))
    do n=1,ntracetloc
#else
    offsetn = 0
    do n=1,ntracet
#endif
       do l = 1,lm(region)
          do j = 1,jm(region)
             do i = 1,im(region)
                rm(i,j,l,n) = rm(i,j,l,n)*m(i,j,l)/fscale(offsetn+n)
#ifdef slopes
                rxm(i,j,l,n) = 0.0
                rym(i,j,l,n) = 0.0
                rzm(i,j,l,n) = 0.0
#ifdef secmom
                rxxm(i,j,l,n) = 0.0
                rxym(i,j,l,n) = 0.0
                rxzm(i,j,l,n) = 0.0
                ryym(i,j,l,n) = 0.0
                ryzm(i,j,l,n) = 0.0
                rzzm(i,j,l,n) = 0.0
#endif
#endif
             end do
          end do
       end do
    end do

    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
    !WP! checked rm with the values in savehdf after scatter 
    !    and gather on 4 processors=>OK

  end subroutine readhdfmmr

  subroutine readhdfmmix(region,file_name)
    !
    ! data are expected in mixing ratio from an mmix file
    !
    use dims,        only : im, jm, lm, nregions, datadir, region_name, parent
    use io_hdf
    use tracer_data, only : mass_dat
    use meteodata  , only : m_dat
    use chem_param,  only : fscale, ntracet, names
    use toolbox,     only : escape_tm
#ifdef MPI
    use mpi_comm,    only : scatter_after_read_t
    use mpi_const,   only : ntracet_ar, ntracetloc
    use mpi_const,   only : com_trac, mpi_integer, root_t, ierr
#endif
    use ParTools,    only : root, myid, root_t

    implicit none

    ! in/out
    integer, intent(in)             :: region
    character(len=*), intent(in)    :: file_name 

    ! local
    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,offsetn, from_file, my_parent
    real,dimension(:,:,:,:), allocatable :: field
    real,dimension(:,:,:),   allocatable :: field3d

    ! start

    m => m_dat(region)%data
    rm => mass_dat(region)%rm_t
#ifdef slopes
    rxm => mass_dat(region)%rxm_t
    rym => mass_dat(region)%rym_t
    rzm => mass_dat(region)%rzm_t
#endif
    rm = 0.0
#ifdef slopes
    rxm = 0.0
    rym = 0.0
    rzm = 0.0
#endif
    allocate(field(-1:im(region)+2,-1:jm(region)+2,1:lm(region),1:ntracet))
    field = 0.0    !09/2004 cmk: to initialise the edges!
    if ( myid == root_t ) 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
             call escape_tm('readhdfmmix : Could not open file and no parent '// &
                  'available :'//file_name )
          else
             from_file = 0
             print *, 'readhdf: Trying to initialise '//trim(region_name(region))// &
                  ' from parent....'
          endif
       endif   
       if ( from_file == 1 ) then   ! read from file:
          print*,' '
          print*,'readhdfmmix: ',file_name,'... opened for READ access.' 

       istat = sffinfo(io, n_datasets, n_file_attrs)

       print*, 'readhdfmmix: Number of datasets and attributes',  &
            n_datasets, n_file_attrs

       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
             print *, 'Readhdfmmix:', names(n), '  not found in dataset --> 0.0 '
             field(1:im(region), 1:jm(region), 1:lm(region), n) = 0.0
          endif
       enddo   

       istat = sfend(io)
       if (istat /= FAIL) then
          print*,'readhdfmmix: ... file closed'
          print*,' '
       else
          call escape_tm('readhdfmmix: ERROR in restart from HDF file')
       end if

       deallocate(field3d) ! no longer needed
       endif  ! from file
    endif ! root_t   
             
#ifdef MPI
    ! broadcast from_file !
    call mpi_bcast(from_file, 1, mpi_integer , root_t, com_trac, ierr)
#endif
    if ( from_file == 1 ) then 
#ifdef MPI
      call scatter_after_read_t(field,im(region),jm(region),lm(region),&
           2,2,0,rm,root)
#else
      rm = field
#endif
#ifdef MPI
      offsetn = sum(ntracet_ar(0:myid-1))
      do n=1,ntracetloc
#else
      offsetn = 0
      do n=1,ntracet
#endif
         do l = 1,lm(region)
            do j = 1,jm(region)
               do i = 1,im(region)
                  rm(i,j,l,n) = rm(i,j,l,n)*m(i,j,l)/fscale(offsetn+n)  ! from mmix to kg tracer
#ifdef slopes
                  rxm(i,j,l,n) = 0.0
                  rym(i,j,l,n) = 0.0
                  rzm(i,j,l,n) = 0.0
#endif
               end do
            end do
         end do
      end do

      deallocate(field)
      call init_child_edge(region)   ! get edge info from parent....
    else   ! initialise from parent:
       call init_child(region)   ! get info from parent....
    end if

    nullify(m)
    nullify(rm)
#ifdef slopes
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
#endif
    !WP! checked rm with the values in savehdf after scatter 
    !    and gather on 4 processors=>OK

  end subroutine readhdfmmix



  !------------------------------------------------------------------
  ! save all essential model parameters and fields on unit kdisk
  !
  !                                       v 8.4
  !  modified for HDF output
  !------------------------------------------------------------------

  subroutine SaveHDF( msg, file_name, status )

    use file_hdf   , only : THdfFile, Init, Done, WriteAttribute, TSds, SetDim, WriteData
    use dims       , only : splitsave
    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 : cdebug, kmain, kdebug
    use dims       , only : limits
    use dims       , only : at, bt
    use dims       , only : nsplitsteps, splitorder
    use tracer_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 toolbox,     only : escape_tm
    use partools   , only : myid, root_t, root_k
    use partools   , only : Par_Barrier, Par_Gather_Tracer_t
    use partools   , only : Par_Gather_Tracer_k
    use partools   , only : Par_Barrier
    use partools   , only : which_par, previous_par
    use global_data, only : outdir
    
    ! --- in/out ----------------------------------

    character(len=*), intent(in)  ::  file_name
    character(len=*), intent(in)  ::  msg
    integer, intent(out)          ::  status
    
    ! --- const ----------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/SaveHDF'
    
    ! --- 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(:,:,:), pointer         ::  m
    real, dimension(:,:,:,:), allocatable   ::  fieldglob4d
    
    type(THdfFile)    ::  hdf
    type(TSds)        ::  sds

    integer           :: region, imr,jmr, lmr,i

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

    if ( .not. splitsave ) call escape_tm('savehdf: splitsave should be true ')

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

    do region = 1,nregions

     if ( myid == root_t ) 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:

      call Par_Barrier

      m => m_dat(region)%data
      rm => mass_dat(region)%rm_t
#ifdef slopes
      rxm => mass_dat(region)%rxm_t
      rym => mass_dat(region)%rym_t
      rzm => mass_dat(region)%rzm_t
#ifdef secmom
      rxxm => mass_dat(region)%rxxm_t
      rxym => mass_dat(region)%rxym_t
      rxzm => mass_dat(region)%rxzm_t
      ryym => mass_dat(region)%ryym_t
      ryzm => mass_dat(region)%ryzm_t
      rzzm => mass_dat(region)%rzzm_t
#endif
#endif

      which_par=previous_par(region)
      if ( which_par /= 'tracer' ) &
           call escape_tm('savehdf: Data should be parallel over tracers')

      imr = im(region) ; jmr = jm(region) ; lmr = lm(region)

      ! allocate global array on all PEs
      allocate(fieldglob4d(-1:imr+2,-1:jmr+2,lmr, ntracet))
      fieldglob4d = 0.0

      ! false signals: do not broadcast
      call Par_Gather_Tracer_t(fieldglob4d,imr,jmr,lmr,2,2,0,ntracet,rm,.false.)

      call Par_Barrier

      if ( myid == root_t ) then 

        call Init( sds, hdf, 'm', shape(m(1:imr,1:jmr,1: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,m(1:imr,1:jmr,1:lmr),status)
        call Done(sds,status)
        IF_ERROR_RETURN(status=1)

        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,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status)
        call Done(sds,status)
        IF_ERROR_RETURN(status=1)

      end if

      call Par_Barrier

#ifdef slopes
      ! false signals: do not broadcast
      call Par_Gather_Tracer_t(fieldglob4d,imr,jmr,lmr,2,2,0,ntracet,rxm,.false.)

      if ( myid == root_t ) 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,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status)
        call Done(sds,status)
        IF_ERROR_RETURN(status=1)

      end if

      call Par_Barrier

      ! false signals: do not broadcast
      call Par_Gather_Tracer_t(fieldglob4d,imr,jmr,lmr,2,2,0,ntracet,rym,.false.)

      if ( myid == root_t ) 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,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status)
        call Done(sds,status)
        IF_ERROR_RETURN(status=1)

      end if
      call Par_Barrier
      ! false signals: do not broadcast
      call Par_Gather_Tracer_t(fieldglob4d,imr,jmr,lmr,2,2,0,ntracet,rzm ,.false.)

      if ( myid == root_t ) 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,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status)
        call Done(sds,status)
        IF_ERROR_RETURN(status=1)

      end if

#ifdef secmom
      ! second moments
      call Par_Barrier
      call Par_Gather_Tracer_t(fieldglob4d,imr,jmr, lmr, 2,2,0,ntracet, rxxm ,.false.)   ! false signals: do not broadcast
      if ( myid == root_t ) 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,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status)
           call Done(sds,status)
           IF_ERROR_RETURN(status=1)
      end if

      call Par_Barrier
      call Par_Gather_Tracer_t(fieldglob4d,imr,jmr, lmr, 2,2,0,ntracet, rxym ,.false.)   ! false signals: do not broadcast

      if ( myid == root_t ) 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,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status)
           call Done(sds,status)
           IF_ERROR_RETURN(status=1)
      end if
      call Par_Barrier
      call Par_Gather_Tracer_t(fieldglob4d,imr,jmr, lmr, 2,2,0,ntracet, rxzm ,.false.)   ! false signals: do not broadcast

      if ( myid == root_t ) 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,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status)
           call Done(sds,status)
           IF_ERROR_RETURN(status=1)
      end if

      call Par_Barrier
      call Par_Gather_Tracer_t(fieldglob4d,imr,jmr, lmr, 2,2,0,ntracet, ryym ,.false.)   ! false signals: do not broadcast

      if ( myid == root_t) 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,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status)
           call Done(sds,status)
           IF_ERROR_RETURN(status=1)
      end if

      call Par_Barrier
      call Par_Gather_Tracer_t(fieldglob4d,imr,jmr, lmr, 2,2,0,ntracet, ryzm ,.false.)   ! false signals: do not broadcast

      if ( myid == root_t) 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,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status)
           call Done(sds,status)
           IF_ERROR_RETURN(status=1)
      endif
      call Par_Barrier
      call Par_Gather_Tracer_t(fieldglob4d,imr,jmr, lmr, 2,2,0,ntracet, rzzm ,.false.)   ! false signals: do not broadcast

      if ( myid == root_t) 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,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status)
           call Done(sds,status)
           IF_ERROR_RETURN(status=1)
       endif

#endif
#endif

      ! clear
      deallocate(fieldglob4d)
      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

      ! 
      ! short lived tracers
      !

      ! only if there are non-transported tracers
      if ( ntrace_chem > 0 ) then

        ! ahoc ...
        if ( root_t /= root_k ) then
          write (gol,'("roots for t and k not equal: ",2i4)') root_t, root_k; call goErr
          write (gol,'("in ",a)') rname; status=1; call goErr
        endif

        ! allocate global array on all PEs
        allocate(fieldglob4d(imr,jmr,lmr,ntrace_chem))

        ! tracers parallel over levels:
        rmc => chem_dat(region)%rm_k

        ! collect on root_k
        call Par_Gather_Tracer_k(fieldglob4d,imr,jmr,lmr,0,0,0,ntrace_chem,rmc,.false.,status)
        IF_NOTOK_RETURN(status=1)

        if ( myid == root_k ) then

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

        end if
        call Par_Barrier

        ! clear
        nullify(rmc)      
        deallocate(fieldglob4d)

      end if  ! ntrace_chem > 0

      !
      ! close
      !

      if ( myid == root_t ) then
        call Done( hdf, status )
        IF_NOTOK_RETURN(status=1)
      end if
      call Par_Barrier

    end do  !regions...

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


    call tstamp(kmain,itau,msg)

    if (cdebug) then
       call tstamp(kdebug,itau,'savehdf')
    end if
    
    ! ok
    status = 0

  end subroutine SaveHDF



end module io_save
