#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_NOTOK_MDF(action) if (status/=0) then; TRACEBACK; action; call MDF_CLose(fid,status); status=1; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
#include "tm5.inc"
!----------------------------------------------------------------------------
!                  TM5                                                      !
!----------------------------------------------------------------------------
!BOP
!
! !MODULE:  Restart_remap
!
! !DESCRIPTION: Read 'rm' from restart files at one vertical resolution,
!               use fill3d to expand to higher (current) vertical resolution.
!\\
!\\
! !INTERFACE: 
!
module Restart_remap
  !
  ! !USES:
  !  
  use GO     , only : gol, goPr, goErr
  use dims   , only : nregions
  
  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public  ::  Restart_remap_rm
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'Restart_remap'
  integer                      ::  fid ! file id for IF_NOTOK_MDF macro
  !
  ! !REVISION HISTORY: 
  !    3 Jan 2012 - Andy Jacobson, borrowed wholesale from Le Sager restart.F90
  !                 most comments and code are his.
  !
  ! !REMARKS:
  ! (1) when an error occurs when accessing MDF files, you should first close
  !   the file before returning. The IF_NOTOK_MDF macro takes care of that.
  !   The only thing you need is to call it like that :
  !   
  !      IF_NOTOK_MDF(fid=xxxx)
  !
  !   where you replace xxxx with the integer id (file handler) of the file
  !   you are accessing. 
  ! 
  !EOP
  !------------------------------------------------------------------------

contains


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


  subroutine Restart_remap_rm( status, region )

    use GO         , only : TrcFile, Init, Done, ReadRc
    use global_data, only : rcfile
    use global_data, only : outdir
    use meteodata  , only : lli
    use GO          , only : Set
    use GO          , only : goMatchValue
    use dims        , only : nregions, at, bt
    use dims        , only : istart, idate, idatei
    !    use dims        , only : nlon360, nlat180, iglbsfc
    use chem_param  , only : ntracet, ntrace_chem, ntrace
    use chem_param  , only : names, tracer_name_len
    use global_data , only : rcfile
    use partools    , only : localComm
    use partools    , only : myid, root
    use partools    , only : ntracetloc, offsetn, tracer_active, tracer_loc
    use partools    , only : lmloc, offsetl
    use partools    , only : MPI_INFO_NULL
    use partools    , only : previous_par
    use tracer_data , only : mass_dat, chem_dat
    use meteodata   , only : lli, levi
    use meteodata   , only : sp_dat, phlb_dat, m_dat
    !use meteodata   , only : slhf_dat, sshf_dat
    use MDF         , only : MDF_Open,           MDF_Close,            MDF_Inquire_Dimension
    use MDF         , only : MDF_Inq_VarID,      MDF_Inquire_Variable, MDF_Inq_DimID
    use MDf         , only : MDF_Var_Par_Access, MDF_INDEPENDENT
    use MDF         , only : MDF_Get_Att
    use MDF         , only : MDF_Get_Var
    use MDF         , only : MDF_READ, MDF_NETCDF4
    ! for restart_remap (ARJ Jan 2012)
    !    use meteo, only : pressure_to_mass
    use Grid        , only : TllGridInfo, TLevelInfo, Init, Done
    use Grid        , only : HPressure
    use Grid        , only : Fill3D

    !
    ! !REVISION HISTORY: 
    !    3 Jan 2012 - Andy Jacobson - remap rm
    !    8 Apr 2011 - P. Le Sager - use IF_NOTOK_MDF macro
    !   28 Apr 2011 - P. Le Sager - Check on tracer availability in restart file.
    !                             - Allows for more tracers in restart file than needed
    !   10 May 2011 - P. Le Sager - Added deallocate statement to work with zoom regions
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    character(len=*), parameter                   ::  rname = 'Restart_remap_rm'
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)                          ::  status
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in), optional                 ::  region
    ! --- local --------------------------------
    type(TrcFile)                                 ::  rcF    
    character(len=256)                            ::  rs_read_dir
    logical                                       ::  exist
    integer                                       ::  imr, jmr, lmr
    integer                                       ::  n, region1, region2
    character(len=256)                            ::  fname
    integer                                       ::  ncid
    integer                                       ::  varid_names
    integer                                       ::  iname
    character(len=tracer_name_len), allocatable   ::  values_names(:)
    integer                                       ::  itr, itr_loc, itr_file
    integer                                       ::  varid_rm
    integer                                       ::  ntracet_restart, dimid
    integer                                       ::  shp(2)
    integer                                       ::  lmr_restart, varid_at, varid_bt
    real, allocatable                             ::  at_restart(:),bt_restart(:)
    type(TllGridInfo)                             ::  lli_restart
    type(TLevelInfo)                              ::  levi_restart
    real, allocatable                             ::  inbuf(:,:,:)

    ! --- begin --------------------------------
    write (gol,'("[restart_remap_rm] ")'); call goPr
    !
    ! read settings from rcfile
    !

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

    ! correct istart ?
    if ( istart /= 61 ) then
       write (gol,'("  - skip; istart not 61 but ",i3)') istart; call goPr
       status=0; return
    end if

    ! init time ?
    if ( any( idate /= idatei ) ) then
       write (gol,'("  - skip; idate not idatei but ",i4,5i2.2)') idate; call goPr
       status=0; return
    end if

    ! region range:
    if ( present(region) ) then
       region1 = region
       region2 = region
    else
       region1 = 1
       region2 = nregions
    end if

    ! loop over ns
    do n = region1, region2

       ! grid size
       imr = lli(n)%nlon
       jmr = lli(n)%nlat
       lmr = levi%nlev

       ! name of restart file
       call Restart_FileName( n, fname, status, dir=trim(rs_read_dir) )
       IF_NOTOK_RETURN(status=1)

       write (gol,'("  - reading tracer mass from ",a)') trim(fname); call goPr

       ! test ...
       inquire( file=fname, exist=exist )
       if ( .not. exist ) then
          write (gol,'("  - restart file not found : ",a)') trim(fname); call goErr
          TRACEBACK; status=1; return
       end if

       ! o open netcdf file

#ifdef MPI      
       ! open existing file:
       call MDF_Open( trim(fname), MDF_NETCDF4, MDF_READ, ncid, status, &
            mpi_comm=localComm, mpi_info=MPI_INFO_NULL )
       IF_NOTOK_RETURN(status=1)
#else
       ! open existing file:
       call MDF_Open( trim(fname), MDF_NETCDF4, MDF_READ, ncid, status )
       IF_NOTOK_RETURN(status=1)
#endif

       ! check for vertical dimension compatibility
       call MDF_Inq_DimID( ncid, 'lev', dimid, status )
       IF_NOTOK_MDF(fid=ncid)

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

       write (gol,'("  - [restart_remap] Model compiled for ",i3," levels; restart file has ",i3" levels.")') lmr, lmr_restart; call goPr

       ! init levi_restart
       allocate(at_restart(lmr_restart+1))
       allocate(bt_restart(lmr_restart+1))
       !
       call MDF_Inq_VarID( ncid, 'at', varid_at, status )
       IF_NOTOK_MDF(fid=ncid)
       !
       call MDF_Get_Var( ncid, varid_at, at_restart(1:(lmr_restart+1)), status )
       IF_NOTOK_MDF(fid=ncid)
       !
       call MDF_Inq_VarID( ncid, 'bt', varid_bt, status )
       IF_NOTOK_MDF(fid=ncid)
       !
       call MDF_Get_Var( ncid, varid_bt, bt_restart(1:(lmr_restart+1)), status )
       IF_NOTOK_MDF(fid=ncid)
       !
       call Init( levi_restart, lmr_restart, at_restart, bt_restart, status )
       IF_NOTOK_RETURN(status=1)

       ! init lli_restart
       call Init( lli_restart, lli(n)%lon_deg(1), lli(n)%dlon_deg, lli(n)%im, &
            lli(n)%lat_deg(1), lli(n)%dlat_deg, lli(n)%jm, status )
       IF_NOTOK_RETURN(status=1)

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

       ! get dimension of "names"
       call MDF_Inquire_Variable( ncid, varid_names, status, shp=shp )
       IF_NOTOK_MDF(fid=ncid)

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

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

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

       ! set independent data mode:
       call MDF_Var_Par_Access( ncid, varid_rm , MDF_INDEPENDENT, status )
       IF_NOTOK_MDF(fid=ncid)

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

       allocate(inbuf(imr,jmr,lmr_restart))

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

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

          if ( status < 0 ) then
             write(gol,'("Requested tracer `",a,"` not found in restart file!")') trim(names(itr)); call goErr
             IF_NOTOK_MDF(fid=ncid)
          endif

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

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

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

             case ( 'tracer' )

                ! on this pe ?
                if ( tracer_active(itr) ) then

                   ! local tracer index:
                   itr_loc = tracer_loc(itr)


                   call MDF_Get_Var( ncid, varid_rm, inbuf, &
                        start=(/1,1,1,itr_file/), &
                        count=(/imr,jmr,lmr_restart,1/), status=status )
                   IF_NOTOK_MDF(fid=ncid)

                   call Fill3D( lli(n), levi, 'n', sp_dat(n)%data(1:imr,1:jmr,1), &
                        mass_dat(n)%rm_t(1:imr,1:jmr,1:lmr,itr_loc), &
                        lli_restart, levi_restart, inbuf(1:imr,1:jmr,1:lmr_restart), 'sum', status )
                   IF_NOTOK_RETURN(status=1)

#ifdef slopes
                   ! Ensure that slopes are initialized to "unset" values of 0.0.
                   ! Wouter says that we could remap levels for rxm et al., but
                   ! 0.0 will also work.  The noise induced from remapping the
                   ! rm array is almost certainly bigger than any issues from
                   ! having this "default=0.0" slopes information.  -ARJ 1 Jan 12

                   mass_dat(n)%rxm_t(1:imr,1:jmr,1:lmr,itr_loc:itr_loc) = 0.0                     
                   mass_dat(n)%rym_t(1:imr,1:jmr,1:lmr,itr_loc:itr_loc) = 0.0                     
                   mass_dat(n)%rzm_t(1:imr,1:jmr,1:lmr,itr_loc:itr_loc) = 0.0                     
#endif

                end if   ! tracer active

             case default

                write (gol,'("remapping only supported for distribution of tracer by PE (not levels or otherwise).")'); call goPr
                write (gol,'("in ",a)') rname; call goErr; status=1; return

             end select
          end if
       end do  ! tracers

       ! free mem for next region
       DeAllocate( values_names)      
       deallocate(inbuf)
       deallocate(at_restart)
       deallocate(bt_restart)


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

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

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

    end do   ! regions

    ! ok
    status = 0

  end subroutine Restart_remap_rm


  ! ***
  ! direct copy of subroutine from restart module. - ARJ
  subroutine Restart_FileName( region, fname, status, key, dir, isfirst )

    use dims       , only : idate
    use global_data, only : outdir
    use meteodata  , only : lli

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

    integer, intent(in)             ::  region
    character(len=*), intent(out)   ::  fname
    integer, intent(out)            ::  status

    character(len=*), intent(in), optional  ::  dir    
    character(len=*), intent(in), optional  ::  key
    logical, intent(in), optional           ::  isfirst

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

    character(len=*), parameter  ::  rname = 'Restart_FileName'
    character(len=256)           ::  adir
    character(len=32)            ::  akey

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

    ! destination directory:
    adir = trim(outdir)
    if ( present(dir) ) adir = trim(dir)

    ! extra key, for example '_x' to denote that 
    ! a restart file was dumped after process 'x':
    akey = ''
    if ( present(key) )  akey = trim(key)

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

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

    ! ok
    status = 0

  end subroutine Restart_FileName

end module Restart_remap
