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

  use GO, only : gol, goErr, goPr
  
  implicit none

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

  private

  public :: default, exitus, start


  ! --- const --------------------------------
    
  character(len=*), parameter  ::  mname = 'initexit'

  
contains


  subroutine default
    !---------------------------------------------------------------------
    ! set default values of control variables
    !
    !                             mh, 27-jun-1989
    !                             mh, 26-sep-1992
    !                             v 9.1
    !                              mk, 21-dec-2002
    !---------------------------------------------------------------------
    use dims
    use datetime,    only : tau2date
    use toolbox,     only : escape_tm
    implicit none

    ! local
    integer :: i,k,n

    ! start
    !                 ! set default calendar 
    !                 ! (real calendar, iyear0=1980)
    icalendo=2

    !                 ! default time steps of basic tasks
    nstep=0
    ndyn=3600*3
    nconv=3600*3
    ndiff=24*3600*1000 !never happes in one month
    ntrans=0
    ndiag=4*3600
    nchem=0
    nsrce=24*3600
    !nread=3*3600       ! <--- set in main program by Meteo_Init
    nwrite=-1
    ninst=0
    !c                       ! default is restart
    istart=10
    itaui=0
    newsrun=.true.
    call tau2date(itaui,idatei)
    if ( mod(idatei(4),3) /= 0 ) &
         call escape_tm('default: GMT start time should be multiple of 3')    
    itaue=itaui
    call tau2date(itaue,idatee)
    itaut=0
    call tau2date(itaut,idatet)


    !c                       ! output for conservation diagnostics
    !c                       ! -1: daily, -2: monthly, -3: yearly,
    !c                               ! >=0 interval in sec
    ndiagp1=-2

    !c                       ! output for mean field diagnostics
    !c                       ! -1: daily, -2: monthly, -3: yearly,
    !c                               ! >=0 interval in sec
    ndiagp2=-2

    !c                               ! full convection
    czeta=1.
    !c                               ! full vertical diffusion
    czetak=1.
    !c                               ! scaling factor for horizontal diffusion
    limits=.true.

    !c                       ! checking interval
    ncheck=6

    !c                       ! control for debug output
    cdebug=.false.
    okdebug=.false.
    revert = 1

  end subroutine default



  !--------------------------------------------------------------
  ! terminate a model run
  !
  ! icode=0   normal termination
  !       else      abnormal termination
  !
  ! msg       character string to be printed on unit kmain
  !
  !                             mh, 29-jun-1989
  !
  ! changed to v 9.1                        mh, 26-sep-1992
  !
  !--------------------------------------------------------------

  subroutine exitus

    use GO           , only : gol, goErr, goPr
    use GO           , only : goGetFU
    use dims         , only : nregions
    use dims         , only : nstep, itau
    use dims         , only : cpu0
    use dims         , only : nxi, xi, nloop_max
    use dims         , only : okdebug, kmain
    use chem_param   , only : ntrace_chem
    use tracer_data  , only : chem_dat
    use global_data  , only : free_fields
    use io_save,       only : savehdf
    use restart,       only : rs_write
    use user_output,   only : User_output_Done
    !use user_output,   only : user_output_done
    use sources_sinks, only : trace_end
    use datetime,      only : tstamp
    use ParTools     , only : myid, root, Par_Barrier, Par_StopMPI
    use tracer_data  , only : Par_Check_Domain
    use global_data  , only : outdir
    use advectm_cfl,    only : done_cfl
#ifdef with_budgets
    use budget_global, only : done_budget_global
#endif

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

    character(len=*), parameter ::  rname = mname//'/exitus'
    
    ! --- local ------------------------------
    
    integer   ::  status
    integer   ::  fu
    integer   ::  region, gatheri
    real      ::  cpu3, gatherr

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

    do region=1,nregions      
       ! transfer data to tracer domain
       call Par_Check_Domain( region, 'n', 'tracer' ) 
    end do

#ifndef with_prism
    ! save the model state.
    ! in case of coupled system, do not use this, as 
    ! restart-files are used. This routine is quite memory-consuming,
    ! and should not be called if not used.
    if(.not.rs_write) then 
       call SaveHDF( 'successful completion of run',trim(outdir)//'/save', status )
       if (status/=0) then; write (gol,'("in ",a,"; continue")') rname; call goErr; end if
    else
       write (gol,'("in ",a,"; no writing of save.hdf files since restart files are used")') rname; call goPr
       write (gol,'("in ",a,"; ....(restart.write : T or with_prism in rc-file)")') rname; call goPr
    endif
#endif

    !>>> moved to TM5/TM5_Model_Run
    !! complete user-specified output   !CMK only once!
    !call user_output_done( status )
    !if (status/=0) then; write (gol,'("in ",a,"; continue")') rname; call goErr; end if
    !<<<

#ifdef with_budgets
    ! save budgets, print sumary
    call done_budget_global ( status )
    if (status/=0) then; write (gol,'("in ",a,"; continue")') rname; call goErr; end if
#endif

    ! finish up sources/sinks:
    call trace_end( status )
    if (status/=0) then; write (gol,'("in ",a,"; continue")') rname; call goErr; end if
    
    ! clear short-lived tracers:
    if ( ntrace_chem > 0 ) then
      do region = 1, nregions
        deallocate( chem_dat(region)%rm_k )
      end do
    end if    

    call Par_Barrier
    if ( okdebug ) then
      write (gol,'("exitus: starting to free fields ",i4)') myid; call goPr
    end if

    call free_fields     !free memory of this run 

    call Par_Barrier
    if ( okdebug ) then
      write (gol,'("exitus: starting time computation",i4)') myid; call goPr
    end if
    call Par_Barrier

    if ( myid == root ) then
      write (gol,'(" ")'); call goPr
      write (gol,'("exitus: CFL info:")'); call goPr
      do region=1, nregions
        write (gol,'(a,3i4,f10.4)') ' exitus: x: region, nxi, nloop_max, xi',  &
               region, nxi(region,1), nloop_max(region,1), xi(region,1); call goPr
        write (gol,'(a,3i4,f10.4)') ' exitus: y: region, nxi, nloop_max, xi',  &
               region, nxi(region,2), nloop_max(region,2), xi(region,2); call goPr
        write (gol,'(a,3i4,f10.4)') ' exitus: z: region, nxi, nloop_max, xi',  &
               region, nxi(region,3), nloop_max(region,3), xi(region,3); call goPr
      end do
    end if
    
    ! cfl finished ...
    call Done_CFL

    call Par_Barrier

    if ( myid == root ) then

      write (gol,'(1x)'); call goPr
      write (gol,'(1x)'); call goPr
      write (gol,'(1x,39("--"))'); call goPr
      write (gol,'(1x,39("--"))'); call goPr
      write (gol,'(1x)'); call goPr
      write (gol,'("exitus: program has terminated normally.")'); call goPr

      call tstamp(kmain,itau,'exitus: final time')

      call cputim(cpu3)
      cpu3 = cpu3-cpu0

      write (gol,'(a,i8,a,f16.2,a,f16.8)') &
           ' exitus: no of timesteps: ',nstep,' cpu (s):',cpu3, &
           ' cpu/step: ',cpu3/nstep; call goPr

      write (gol,'(1x)'); call goPr
      write (gol,'(1x,39("--"))'); call goPr
      write (gol,'(1x,39("--"))'); call goPr
      write (gol,'(1x)'); call goPr

    end if ! all Pe's from here 

    call Par_barrier
    if ( okdebug ) then
      write (gol,'("exitus: calling stopmpi ",i4)') myid; call goPr
    end if
    call Par_barrier
    
    ! end program:
    !  o mpi_finalize
    !  o stop
    !call Par_StopMPI 

  end subroutine exitus
      
      

  !-----------------------------------------------------------------
  !     this subroutine performs the initialization of a model run or 
  !     its continuation.
  !
  !
  !     structure of control input :
  !     line 1 :     blank or name of alternative control input file
  !     line 2-3:    text for run header/info (2*80 char)
  !     line 4-... : namelist inputz
  !                  (must be terminated by string
  !                  '&end' in column 2-5)
  !     subsequent input may be read in sr trace0 (via unit kinput0)
  !
  !                                       mh,  1-dec-1990
  !
  !      changed to v 8.5                 mh, mpi HH     23-feb-1994
  !      changed to zoom version          MK
  !-----------------------------------------------------------------------

  subroutine start( tread1, tread2, status )

    use GO           , only : gol, goPr, goErr, goLabel
    use GO           , only : TrcFile, Init, Done, ReadRc
    use GO           , only : TDate, NewDate, IncrDate, AnyDate
    use GO           , only : operator(+), operator(-), rTotal
    use GO           , only : goTranslate
    use dims         , only : im, jm, lm, dxy11, nlat180, nlon360
    use dims         , only : nregions, region_name
    use dims         , only : ndyn, nchem, nsrce, nconv, nread, ndyn_max
    use dims         , only : ndiag, ndiagp1, ndiagp2, nwrite, ninst, ndiff
    use dims         , only : nstep0, nstep
    use dims         , only : idatei, idatee, idatet, idate, revert
    use dims         , only : itaui, itaue, itaut, itau
    use dims         , only : itaur
    use dims         , only : iyear0, julday0, icalendo
    use dims         , only : newyr, newday, newmonth, newsrun
    use dims,          only : newmonth, idate
    use dims,          only : mlen, sec_day,sec_month,sec_year
    use dims         , only : okdebug, cpu0, kmain
    use dims         , only : istart
    use dims         , only : czeta, czetak, limits, cdebug
    use chem_param   , only : ntracet, ntrace_chem
#ifdef with_pycasso
    use global_data  , only : rcfile
#else
    use global_data  , only : rcfile, rcfile_runtime
#endif
    use global_data  , only : inputdir, outdir
    use global_data  , only : region_dat
    use tracer_data  , only : mass_dat, chem_dat
    use global_data  , only : declare_fields
    use global_data  , only : fcmode, tfcday0
    use Meteo        , only : Meteo_Setup_Mass, Meteo_Setup_Other
    use sources_sinks, only : trace0, trace1
#ifdef with_budgets
    use budget_global, only : Init_budget_global
#endif
    use redgridZoom  , only : grid_reduced, InitRedGrid
    use redgridZoom  , only : nred, nredmax, clustsize, imred, imredj, jred
    use geometry     , only : geomtryh, calc_dxy
    use io_save      , only : readhdf, readhdfmmr, savehdf, readhdfmmix
    use user_output  , only : User_output_Init, user_output_step
    use datetime,      only : calc_sm
    use datetime     , only : date2tau, tau2date, julday, tstamp
    use zoom_tools   , only : coarsen_region, update_parent
    use advect_tools , only : advect_m, m2phlb, m2phlb1, coarsen_ambm
#ifdef MPI    
    use mpi_const    , only : localComm
    use mpi_const    , only : MPI_INTEGER, MPI_LOGICAL, MPI_CHARACTER
    use mpi_const    , only : MY_REAL
    use mpi_const    , only : ierr
#endif
    use ParTools     , only : myid, root, Par_Barrier
    use ParTools     , only : lmloc
    use advectm_cfl  , only : Init_CFL
#ifndef without_advection
    use advectm_cfl  , only : Check_CFL
#endif
    use user_input   , only : user_input_start
#ifdef oasis3
    use tm5_prism    , only : PRISM_start_date
#endif
    use restart      , only : Restart_Read, Restart_Write
    use restart_remap, only : Restart_Remap_rm
#ifdef with_tendencies
#ifdef oasis4
    use tm5_tendency_eval    , only : prism_init_time
#endif
    use tm5_tendency_eval    , only : set_init_tend
#endif

    use meteodata, only : sp_dat
    
    ! --- in/out --------------------------------------------
    
    type(TDate), intent(out)    ::  tread1, tread2
    integer, intent(out)        ::  status

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

    character(len=*), parameter ::  rname = mname//'/start'
    
    ! --- local ---------------------------------------------
    
    integer             ::  n, region
    integer             ::  imr, jmr, lmr
    character(len=256)  ::  fname, fdir
    character(len=32 )  ::  key
    type(TrcFile)       ::  rcF
    type(TDate)         ::  tdyn, tr(2)
#ifdef with_pycasso
    character(len=32)   ::  stime
#else
    character(len=8)    ::  ccyymmdd
#endif
    integer             ::  ccyy, mm, dd

    ! --- begin -----------------------------------------------
    
    call goLabel( rname )

    write (gol,'(a,": begin")') rname; call goPr

    write (gol,'(a,":   default ...")') rname; call goPr

    ! set default control parameters
    call default

    write (gol,'(a,":   rcfile settings ...")') rname; call goPr

    if ( myid == root ) then

      !
      ! read parameters from rcfile
      !

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

      ! grid redudction ?
      call ReadRc( rcF, 'proces.advection.reduced', grid_reduced, status )
      IF_NOTOK_RETURN(status=1)

      ! forecast series ?
      call ReadRc( rcF, 'time.fc', fcmode, status, default=.false. )
      IF_ERROR_RETURN(status=1)
      ! read forecast day 0 ?
      if ( fcmode ) then
#ifdef with_pycasso
        ! read day: yyyy-mm-dd
        call ReadRc( rcF, 'time.fc.day0' , stime, status )
        IF_NOTOK_RETURN(status=1)
        call goTranslate( stime, '/-', ' ', status )
        IF_NOTOK_RETURN(status=1)
        read (stime,*,iostat=status) ccyy, mm, dd
        if ( status /= 0 ) then
          write (gol,'("reading ccyy mm dd from : ",a)') trim(stime); call goErr
          TRACEBACK; call goErr; status=1; return
        end if
#else
        ! read day:
        call ReadRc( rcF, 'time.fc.day0', ccyymmdd, status )
        IF_NOTOK_RETURN(status=1)
        ! extract year, month, day:
        read (ccyymmdd,'(i4.4,2i2.2)',iostat=status) ccyy, mm, dd
        if ( status /= 0 ) then
          write (gol,'("reading ccyy mm dd from : ",a)') ccyymmdd; call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
        end if
#endif
        ! store day:
        tfcday0 = NewDate(year=ccyy,month=mm,day=dd)
      else
        ! dummy:
        tfcday0 = AnyDate()
      end if

      call ReadRc(rcF, 'time.ndyn_max', ndyn, status )
      IF_NOTOK_RETURN(status=1)
      nconv = ndyn
      nchem = ndyn
      nsrce = ndyn
      ndyn_max = ndyn

      ! ensure that every 'nread' seconds is at the end of a dynamic time step:
      call ReadRc( rcF, 'time.ntimestep', nread, status )
      IF_NOTOK_RETURN(status=1)

      ! how to initialize the tracer fields ?
      call ReadRc( rcF, 'istart', istart, status )
      IF_NOTOK_RETURN(status=1)

      ! input files:
      call ReadRc( rcF, 'inputdir', inputdir, status )
      IF_NOTOK_RETURN(status=1)

      ! print debug info ?
      call ReadRc( rcF, 'okdebug', okdebug, status )
      IF_NOTOK_RETURN(status=1)

#ifdef with_pycasso
      ! name of output directory:
      call ReadRc( rcF, 'output.dir', outdir, status )
      IF_NOTOK_RETURN(status=1)
#else
      ! default name for output directory;
      ! might be overwritten later on by name including time range:
      call ReadRc( rcF, 'outputdir', outdir, status )
      IF_NOTOK_RETURN(status=1)
#endif

#ifdef with_pycasso
#else
      ! close
      call Done( rcF, status )
      IF_NOTOK_RETURN(status=1)
#endif

!#ifdef with_prism 
#ifdef oasis4
      call prism_init_time
#else
#ifdef with_pycasso
      ! start time:
      call ReadRc( rcF, 'jobstep.timerange.start' , stime, status )
      IF_NOTOK_RETURN(status=1)
      call goTranslate( stime, '/-:', ' ', status )
      IF_NOTOK_RETURN(status=1)
      read (stime,*,iostat=status) idatei
      if ( status /= 0 ) then
        write (gol,'("could not read start time from : ",a)') trim(stime); call goErr
        TRACEBACK; status=1; return
      end if

      ! end time:
      call ReadRc( rcF, 'jobstep.timerange.end' , stime, status )
      IF_NOTOK_RETURN(status=1)
      call goTranslate( stime, '/-:', ' ', status )
      IF_NOTOK_RETURN(status=1)
      read (stime,*,iostat=status) idatee
      if ( status /= 0 ) then
        write (gol,'("could not read end time from : ",a)') trim(stime); call goErr
        TRACEBACK; status=1; return
      end if

      ! 'target' time?
      idatet = idatee

      ! close:
      call Done( rcF, status )
      IF_NOTOK_RETURN(status=1)
#else
      !
      ! read variables from timing rc file
      !

      call Init( rcF, rcfile_runtime, status )
      IF_NOTOK_RETURN(status=1)

      call ReadRc( rcF, 'year1' , idatei(1), status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'month1', idatei(2), status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'day1'  , idatei(3), status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'hour1' , idatei(4), status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'minu1' , idatei(5), status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'sec1'  , idatei(6), status )
      IF_NOTOK_RETURN(status=1)

      call ReadRc( rcF, 'year2' , idatee(1), status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'month2', idatee(2), status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'day2'  , idatee(3), status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'hour2' , idatee(4), status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'minu2' , idatee(5), status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'sec2'  , idatee(6), status )
      IF_NOTOK_RETURN(status=1)


      idatet = idatee


      call ReadRc( rcF, 'istart', istart, status )
      IF_NOTOK_RETURN(status=1)

      call ReadRc( rcF, 'outputdir', outdir, status )
      IF_NOTOK_RETURN(status=1)

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

#endif
#endif

    end if

    write (gol,'(a,":   broadcast ...")') rname; call goPr

#ifdef MPI      
    ! broadcast namelist
    call MPI_BCAST(istart   ,1,   MPI_INTEGER,    root ,localComm,ierr)  
    call MPI_BCAST(ndyn     ,1,   MPI_INTEGER,    root ,localComm,ierr) 
    call MPI_BCAST(ndyn_max ,1,   MPI_INTEGER,    root ,localComm,ierr)   
    call MPI_BCAST(nconv    ,1,   MPI_INTEGER,    root ,localComm,ierr)  
    call MPI_BCAST(ndiag    ,1,   MPI_INTEGER,    root ,localComm,ierr)  
    call MPI_BCAST(nchem    ,1,   MPI_INTEGER,    root ,localComm,ierr)  
    call MPI_BCAST(nsrce    ,1,   MPI_INTEGER,    root ,localComm,ierr)  
    call MPI_BCAST(nread    ,1,   MPI_INTEGER,    root ,localComm,ierr)  
    call MPI_BCAST(nwrite   ,1,   MPI_INTEGER,    root ,localComm,ierr)  
    call MPI_BCAST(ninst    ,1,   MPI_INTEGER,    root ,localComm,ierr)  
    call MPI_BCAST(ndiff    ,1,   MPI_INTEGER,    root ,localComm,ierr)  
    call MPI_BCAST(icalendo ,1,   MPI_INTEGER,     root ,localComm,ierr)  
    call MPI_BCAST(iyear0   ,1,   MPI_INTEGER,     root ,localComm,ierr)  
    call MPI_BCAST(idatei   ,6,   MPI_INTEGER,     root ,localComm,ierr)  
    call MPI_BCAST(idatee   ,6,   MPI_INTEGER,     root ,localComm,ierr)  
    call MPI_BCAST(idatet   ,6,   MPI_INTEGER,     root ,localComm,ierr)  
    call MPI_BCAST(ndiagp1  ,1,   MPI_INTEGER,     root ,localComm,ierr)  
    call MPI_BCAST(ndiagp2  ,1,   MPI_INTEGER,     root ,localComm,ierr)  
    call MPI_BCAST(czeta    ,1,   MY_REAL,         root ,localComm,ierr)
    call MPI_BCAST(czetak   ,1,   MY_REAL,         root ,localComm,ierr)
    call MPI_BCAST(limits   ,1,   MPI_LOGICAL,     root ,localComm,ierr)
    call MPI_BCAST(cdebug   ,1,   MPI_LOGICAL,     root ,localComm,ierr)
    call MPI_BCAST(okdebug  ,1,   MPI_LOGICAL,     root ,localComm,ierr)
    call MPI_BCAST(    inputdir, len(inputdir), MPI_CHARACTER, root, localComm, ierr )
    call MPI_BCAST(      outdir, len(  outdir), MPI_CHARACTER, root, localComm, ierr )
    call MPI_BCAST(grid_reduced, 1            , MPI_LOGICAL  , root ,localComm, ierr )
#endif

    write (gol,'(a,":   times ...")') rname; call goPr

    !c                             !--------------------------!
    !c                             ! initialize new model run !
    !c                             !--------------------------!
    !c                             ! initialize time parameters
    !c                             ! and set calendar
    !c                             !--------------------------!
    !
    ! itau runs from beginning of year -1 allows running from 1-1-yyyy 
    ! cmk aug/2003
    iyear0 = idatei(1)-1 
    julday0=julday(1,1,iyear0) 
    if (icalendo.eq.2) then
       julday0=julday(1,1,iyear0)
    end if

    call date2tau(idatei,itaui)
    itau=itaui
    call tau2date(itau,idate)
    call declare_fields           !allocate memory for this run

    call date2tau(idatee,itaue)
    call date2tau(idatet,itaut)

    ! set time flags
    newyr=.true.
    newday=.true.
    newmonth=.true.
    newsrun=.true.

    ! step counter
    nstep0=0
    nstep=0

#ifdef with_prism
#ifdef oasis3
      ! store initial time for prism coupling via oasis3
      ! (time is defined as seconds from begin)
      PRISM_start_date = idatei
      !
#endif
#endif

    write (gol,'(a,":   init cfl ...")') rname; call goPr

    ! initialise CFL:
    CALL init_cfl          !initialise CFL preventing routine

    write (gol,'(a,":   determine children etc ...")') rname; call goPr

    ! cmk: DO THIS FIRST: CHILDREN ARE NEEDED FOR REDGRID TEST!
    call determine_children_etc( status )
    IF_NOTOK_RETURN(status=1)

    write (gol,'(a,":   reduced grid ...")') rname; call goPr

    ! set up model geometry...
    if ( grid_reduced ) then
    
      do n = 1, nregions
        if ( myid == root ) then
          call initredgrid( n, im(n), jm(n), status )
          IF_NOTOK_RETURN(status=1)
        end if
#ifdef MPI      
        call MPI_BCAST( status, 1, MPI_INTEGER, root, localComm, ierr )
#endif
        IF_NOTOK_RETURN(status=1)
      end do

#ifdef MPI      
      call MPI_BCAST( nred     , nregions         , MPI_INTEGER, root, localComm, ierr)
      call MPI_BCAST( clustsize, nredmax*nregions, MPI_INTEGER, root, localComm, ierr)
      call MPI_BCAST( imred    , 180             , MPI_INTEGER, root, localComm, ierr)
      call MPI_BCAST( imredj   , nredmax*nregions, MPI_INTEGER, root, localComm, ierr)
      call MPI_BCAST( jred     , nredmax*nregions, MPI_INTEGER, root, localComm, ierr)
#endif

    end if

    write (gol,'(a,":   split order ...")') rname; call goPr

    call define_splitorderzoom

    write (gol,'(a,":   blank zoom regions ...")') rname; call goPr

    do region = 1,nregions
       call blank_zoom_region(region,region_dat(region))
    end do

    !! display zoom regions:
    !if ( myid == root ) call print_zoom(region_dat)

    
    write (gol,'(a,":   horizontal geometry ...")') rname; call goPr

    ! grid cell area in 1x1 grid
    call calc_dxy(dxy11,nlat180)

    ! grid definition
    do region = 1, nregions
       call GeomtryH( region )
    end do
    

    write (gol,'(a,":   initial meteo ...")') rname; call goPr

    ! setup from start time to end of interval  [k*nread,(k+1)*nread]
    tread1 = NewDate(time6=idate)
    tread2 = NewDate(time6=(/idate(1:3),0,0,0/)) + IncrDate(sec=floor(idate(4)*3600.0/nread)*nread+nread)

    ! n is the number of dynamic intervals within the
    ! time interval for which the meteo has been setup:
    n = ceiling( rTotal(tread2-tread1,'sec') / real(ndyn) )
    ndyn = nint( rTotal(tread2-tread1,'sec') / n )

    ! setup pressure and mass fields
    !  o do not check pressure implied by advection
    call Meteo_Setup_Mass( tread1, tread2, status, isfirst=.true. )
    IF_NOTOK_RETURN(status=1)

#ifndef without_advection
    ! determine dynamic timestep ndyn for this interval [tread1,tread2] ;
    ! the initial number of time steps n is increased until no cfl
    ! violations occure
    call Check_CFL( tread1, tread2, n, status )
    IF_NOTOK_RETURN(status=1)
#endif

    ! * setup meteo for dynamic step tdyn+[0,ndyn]
    ! current time (begin of dynamics step)
    tdyn = NewDate( time6=idate  )
    ! time range of dynamic step:
    tr(1) = tdyn
    tr(2) = tdyn + IncrDate( sec=ndyn )
    !! convert pu/pv to am/bm/cm, eventually time interpolated
    !call Setup_MassFlow( tr, status )
    !if (status/=0) call escape_tm('ERROR in tracer') 
    ! setup (interpolate?) other meteo:
    call Meteo_Setup_Other( tr(1), tr(2), status, isfirst=.true. )
    IF_NOTOK_RETURN(status=1)


    !
    ! ** init tracer fiels ****************************
    ! 

    
    write (gol,'(a,":   init tracer fields ...")') rname; call goPr

    ! allocate non transported tracers:
    if ( newsrun ) then
      do region = 1, nregions 
        if ( ntrace_chem > 0 ) then
          imr = im(region) ; jmr = jm(region) ; lmr = lm(region)
          allocate( chem_dat(region)%rm_k(imr,jmr,lmloc,ntracet+1:ntracet+ntrace_chem) )
          chem_dat(region)%rm_k(:,:,:,:) = 0.0
        else
          nullify( chem_dat(region)%rm_k )
        end if
      end do

    end if

    ! get initial tracer fields... (m is also available...)
    IF(revert ==1)THEN

      select case (istart)
        case(1)

          ! initial tracer fields are set
          ! to zero
          do region = 1, nregions
#ifdef MPI          
            mass_dat(region)%rm_k = 0.0
#endif
            mass_dat(region)%rm_t = 0.0

#ifdef slopes
#ifdef MPI           
            mass_dat(region)%rxm_k=0.
            mass_dat(region)%rym_k=0.
            mass_dat(region)%rzm_k=0.
#endif
            mass_dat(region)%rxm_t=0.
            mass_dat(region)%rym_t=0.
            mass_dat(region)%rzm_t=0.
#ifdef secmom
#ifdef MPI           
            mass_dat(region)%rxxm_k=0.
            mass_dat(region)%rxym_k=0.
            mass_dat(region)%rxzm_k=0.
            mass_dat(region)%ryym_k=0.
            mass_dat(region)%ryzm_k=0.
            mass_dat(region)%rzzm_k=0.
#endif
            mass_dat(region)%rxxm_t=0.
            mass_dat(region)%rxym_t=0.
            mass_dat(region)%rxzm_t=0.
            mass_dat(region)%ryym_t=0.
            mass_dat(region)%ryzm_t=0.
            mass_dat(region)%rzzm_t=0.
#endif
#endif
          end do

        !
        ! initial tracer fields are calculated in trace1
        !

        case ( 2 )

          call trace1

          ! transfer information from fine--->coarse...      
          write (gol,'(a,":   update parents ...")') rname; call goPr
          do region = nregions,2,-1
            call update_parent( region )
          end do

        !
        ! initial tracer fields are obtained
        ! from a save file of a previous
        ! model run. m is also obtained!
        !

        case ( 3 )

          !
          ! On istart == 3 (restart) we check if the savefile names were specified
          ! Luckily, this asks for problems, because we also use istart=3 for
          ! rebsubmitting. So we hack in bash scripts with an enormous concentration
          ! of errors. Therefore, I rebuild istart=30 agagin.
          !
           !
           ! Case 3 [pls, 17/9/2010]:
           !
           ! Prior to pycasso : look for files into "savedir" key. This
           ! required a save.cmd and some post-processing.
           !
           
           ! With pycasso: look into "output.dir" (note that save files are
           ! not link/copied/.. to the savedir during post processing anymore:
           ! we leave that kind of work to the user -too many scenarii-).
           ! Corollary: this option should NOT be set by the user, only by a
           ! script when resubmitting a job that do not use restart. And the
           ! savedir key has been made obsolete. Users, who want to start with
           ! save files, must use istart=30 or 31 to specify fully qualified
           ! save filenames. 

#ifdef with_pycasso
          ! open runtime rcfile:
          call Init( rcF, rcfile, status )
          IF_NOTOK_RETURN(status=1)

          ! directory with save files: 
          call ReadRc( rcF, 'prev.output.dir', fdir, status )

          ! close runtime file
          call Done( rcF, status )
          IF_NOTOK_RETURN(status=1)
#else
          ! open rcfile:
          call Init( rcF, rcfile_runtime, status )
          IF_NOTOK_RETURN(status=1)

          ! directory with save files:
          call ReadRc( rcF, 'savedir', fdir, status )

          ! close rcfile:
          call Done( rcF, status )
          IF_NOTOK_RETURN(status=1)
#endif
          
          ! open regular rcfile:
          call Init( rcF, rcfile, status )
          IF_NOTOK_RETURN(status=1)

          ! loop over regions:
          do region = 1, nregions
            ! name of standard save file : save_2000101000_glb3x2.hdf
            write (fname,'(a,"/save_",i4.4,3i2.2,"_",a,".hdf")') &
                trim(fdir), idate(1:4), trim(region_name(region))
            ! read from save file:
            call ReadHDF( region, fname )
          end do

          ! close regular rc file
          call Done( rcF, status )
          IF_NOTOK_RETURN(status=1)

          ! transfer information from fine--->coarse...      
          write (gol,'(a,":   update parents ...")') rname; call goPr
          do region = nregions,2,-1
            call update_parent( region )
          end do

        !
        ! idem, but also able to convert grids and levels
        !
        
        case ( 30 )

           !
           ! On istart == 30, we use save-files defined in the rc-file, with
           ! the following format for the key(s):
           !
           !   start.30.<region-name> : <fully-qualified-filename>
           !
           ! Example:
           ! 
           !   start.30.glb6x4    : /data/save_files/mystuff_glb6x4.hdf
           !

          ! open regular rcfile:
          call Init( rcF, rcfile, status )
          IF_NOTOK_RETURN(status=1)

          ! loop over regions:
          do region = 1, nregions
            ! name of standard save file : save_2000101000_glb3x2.hdf

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

            write (key,'("start.30.",a)') trim(region_name(region))
            call ReadRc( rcF, key, fname, status, default='file_name_empty' )
            IF_NOTOK_RETURN (status = 1)
            write (gol,*) 'Using save file names from rc-file start.30.* values'; call goPr
            ! read from save file:
            call ReadHDF( region, fname )
          end do

          ! close regular rc file
          call Done( rcF, status )
          IF_NOTOK_RETURN(status=1)

          ! transfer information from fine--->coarse...      
          write (gol,'(a,":   update parents ...")') rname; call goPr
          do region = nregions,2,-1
            call update_parent( region )
          end do

        !
        ! idem, but also able to convert grids and levels
        !
          
        case ( 31 )

          call Init( rcF, rcfile, status )
          IF_NOTOK_RETURN(status=1)
          
          do region=1,nregions
          
            ! name of save file for this region:
            call ReadRc( rcF, 'start.31.'//trim(region_name(region)), fname, status )
            IF_NOTOK_RETURN(status=1)
          
            ! read all tracers, convert grid:
            call ReadTracer_saved( region, fname, status )
            IF_NOTOK_RETURN(status=1)
          
            ! overwrite with TM4 fields ?
            call ReadRc( rcF, 'start.31.'//trim(region_name(region))//'.TM4', fname, status, 'none' )
            IF_ERROR_RETURN(status=1)
            ! key found ? then read:
            if ( fname /= 'none' ) then
              call ReadTracer_saved( region, fname, status, tm4=.true. )
              IF_NOTOK_RETURN(status=1)
            end if

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

          ! transfer information from fine--->coarse...      
          write (gol,'(a,":   update parents ...")') rname; call goPr
          do region = nregions,2,-1
            call update_parent( region )
          end do

        !
        ! read from nc restart file:
        !
        
        case ( 33 )
        
          ! read tracer masses and tendencies;
          ! surface pressure is read in meteo
        
          call Restart_Read( status, tracer_mass=.true., tendencies=.true. , air_mass=.true.)
          IF_NOTOK_RETURN(status=1)
          
          ! NOTE: Do not update parents !
          ! A restart file is written after all processes in a time step are finished,
          ! and therefore the parents should be updated alread.
          ! If the 'update_parent' routine is called again here, it causes tiny
          ! differences between a long run and two short ones.

        !
        ! initial tracer fields are obtained
        ! from a save file of a previous
        !

        case ( 4 )

          !mk now HDF
          !WP! read tracer mixing ratios from saveoldfile

          call Init( rcF, rcfile, status )
          IF_NOTOK_RETURN(status=1)
          do region=1,nregions
            call ReadRc( rcF, 'start.4.'//trim(region_name(region)), fname, status )
            IF_NOTOK_RETURN(status=1)
            call readhdfmmr( region, fname )
          end do
          call Done( rcF, status )
          IF_NOTOK_RETURN(status=1)

          ! transfer information from fine--->coarse...      
          write (gol,'(a,":   update parents ...")') rname; call goPr
          do region = nregions,2,-1
            call update_parent( region )
          end do

        !
        ! Initial tracer fields are obtained from a mmix file.
        !
        ! This is typically the choice for combining different versions
        ! or extending the number of tracers
        ! The compounds are searched by name. If not in the mmix file
        ! the field is initialized as zero
        !

        case ( 5 )

          call Init( rcF, rcfile, status )
          IF_NOTOK_RETURN(status=1)
          do region=1,nregions
            call ReadRc( rcF, 'start.5.'//trim(region_name(region)), fname, status )
            IF_NOTOK_RETURN(status=1)
            call readhdfmmix(region,fname)
          end do
          call Done( rcF, status )
          IF_NOTOK_RETURN(status=1)
          
          ! transfer information from fine--->coarse...      
          write (gol,'(a,":   update parents ...")') rname; call goPr
          do region = nregions,2,-1
            call update_parent( region )
          end do

        case(6)   ! start up field from a customized co2 field on 1x1 degree, then coarsened
      
          !call start_co2
          

        ! instart 61:  start up field interpolated from lower vertical resolution fields
        ! saved in a restart file, then interpolated to higher resolution.  m, sp, pblh
        ! are all handled by standard TM5 startup/meteo.F90.

        case( 61 )

          call restart_remap_rm(status)
          IF_NOTOK_RETURN(status=1)
          call restart_write(status=status,key="_remapped")
          IF_NOTOK_RETURN(status=1)

        case(62)   ! start up co2 fields (rm) from restart
      
          write (gol,'(a,":   read tracer masses from restart file ...")') rname; call goPr
          call Restart_Read( status, tracer_mass=.true., tendencies=.false., &
               air_mass=.false., surface_pressure=.false.,pressure=.false., &
               override_istart=.true.)
          IF_NOTOK_RETURN(status=1)
          
        !
        ! User defined
        !
        
        case ( 9 )
        
          call user_input_start( status )
          IF_NOTOK_RETURN(status=1)
        
          ! transfer information from fine--->coarse...      
          write (gol,'(a,":   update parents ...")') rname; call goPr
          do region = nregions,2,-1
            call update_parent( region )
          end do
        
        !
        ! error ...
        !

        case default

          write (gol,'("unsupported istart : ",i2)') istart; call goErr
          call goErr; status=1; return
          
      end select
      
    end if   ! forward run

#ifdef with_tendencies
    call set_init_tend (rcfile,status)
    IF_NOTOK_RETURN(status=1)
#endif

    write (gol,'(a,":   cpu time counting ...")') rname; call goPr

    call cputim(cpu0)
    call Par_Barrier

    write (gol,'(a,":   trace0 ...")') rname; call goPr

    call trace0( status )
    IF_NOTOK_RETURN(status=1)
    call Par_Barrier


#ifdef with_budgets

    write (gol,'(a,":   budgets ...global:")') rname; call goPr
    call init_budget_global ( status )

    IF_NOTOK_RETURN(status=1)

#endif

    !>>> moved to TM5/TM_model_run
    !write (gol,'(a,":   user output ...")') rname; call goPr
    !
    !! initialise user-specified output (!CMK only once!)
    !call user_output_init( status )
    !IF_NOTOK_RETURN(status=1)
    !
    !! write output for start time:
    !itaur = itau   ! <--- seems to be necessary for user_output ...
    !do region = 1, nregions
    !  call user_output_step( region, status )
    !  IF_NOTOK_RETURN(status=1)
    !end do
    !<<<
        
    write (gol,'(a,": model start up complete")') rname; call goPr
    
    ! ok
    call goLabel(); status=0

  end subroutine start


  ! ***

  
  subroutine ReadTracer_saved( region, fname, status, tm4 )

    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 : myid, root
    use partools   , only : previous_par, tracer_active, ntracetloc, tracer_loc
    use partools   , only : Par_Barrier
#ifdef MPI
    use mpi_comm   , only : scatter_after_read_k_3d
#endif
    use tracer_data  , only : mass_dat, chem_dat
    use meteo        , only : lli, levi, sp_dat, m_dat

    ! --- in/out -------------------------------
    
    integer, intent(in)           ::  region
    character(len=*), intent(in)  ::  fname
    integer, intent(out)          ::  status
    
    logical, intent(in), optional     ::  tm4
    
    ! --- const -------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/ReadTracer_saved'
    
    ! 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
    real, dimension(:,:,:,:), pointer   ::  rm

    type(THdfFile)                   ::  hdf
    type(TSds)                       ::  sds_rm, sds_rmc
    integer                          ::  n, n_loc
    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(:,:,:,:)
    real, allocatable                ::  rmt(:,:,:)
    integer                          ::  l
    logical                          ::  with_tm4

    ! --- begin -------------------------------
    
    ! should be parallel over tracers ...
    if ( previous_par(region) /= 'tracer' ) then
      write (gol,'("expected parallel over tracers:")'); call goErr
      write (gol,'("  region       : ",i3)') region; call goErr
      write (gol,'("  previous_par : ",a)') trim(previous_par(region)); call goErr
      TRACEBACK; call goErr; status=1; return
    end if
    
    ! grid sizes
    imr = im(region)
    jmr = jm(region)
    lmr = lm(region)
    
    ! set pointers
    sp =>   sp_dat(region)%data
    m  =>    m_dat(region)%data
    rm => mass_dat(region)%rm_t

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

    ! open hdf file (all pe)
    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) )
    allocate( rmt(imr,jmr,lmr) )
    
    ! 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 .ne. 0) write (gol,'("WARNING - No short lived tracer array found in save file...continueing 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
        write (gol,'("WARNING - tracer not found in save file:")'); call goPr
        write (gol,'("WARNING -   file     : ",a)') trim(fname); call goPr
        write (gol,'("WARNING -   tracer   : ",i3," ",a)') n, trim(names(n)); call goPr
        cycle
      end if

      ! transported or short lived ?
      if ( n <= ntracet ) then
    
        ! not on this proces ?
        if ( .not. tracer_active(n) ) cycle

        ! 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

        ! local index:
        n_loc = tracer_loc(n)

        ! convert grids and levels:
        call Fill3D( lli(region), levi, 'n', sp(1:imr,1:jmr,1), rm(1:imr,1:jmr,1:lmr,n_loc), &
                     lli_in, levi_in, rm_in(:,:,:,1), 'sum', status )
        IF_NOTOK_RETURN(status=1)

        ! info
        write (gol,*) 'initialised tracer ', n, names(n), '; max mass, max vmr : ', &
                           maxval(rm(1:imr,1:jmr,1:lmr,n_loc)), &
                           maxval(rm(1:imr,1:jmr,1:lmr,n_loc)/m(1:imr,1:jmr,1:lmr)*fscale(n)); call goPr
       
      else
      
        ! short lived tracer
        if ( myid == root ) then

          ! 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 

          ! convert grids and levels:
          call Fill3D( lli(region), levi, 'n', sp(1:imr,1:jmr,1), rmt, &
                       lli_in, levi_in, rm_in(:,:,:,1), 'sum', status )
          IF_NOTOK_RETURN(status=1)

          ! info
          write (gol,*) 'initialised tracer ', n, names(n), '; max mass, max vmr : ', &
                           maxval(rmt), &
                           maxval(rmt/m(1:imr,1:jmr,1:lmr)*fscale(n)); call goPr

        else
          ! for safety ..
          rmt = 0.0
        end if
        
        ! distribute over levels:
#ifdef MPI
        call scatter_after_read_k_3d( rmt, imr,jmr,lmr, 0,0,0, &
                                      chem_dat(region)%rm_k(:,:,:,n), root )
#else
        chem_dat(region)%rm_k(:,:,:,n) = rmt
#endif

      end if
      
    end do
    
    ! clear
    deallocate( names_in )
    deallocate( rm_in )
    deallocate( rmt )
    
    ! done with grid definitions:
    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)
    
    ! wait ...
    call Par_Barrier
    
    ! ok
    status = 0
    
  end subroutine ReadTracer_saved
    

  !
  ! ***********************************************************************
  ! cpu time
  ! ***********************************************************************
  !

  subroutine cputim(time )
    implicit none

    real,intent(out) :: time
    
    integer :: clock,clockrate

    call system_clock(clock,clockrate)    
    time=clock*(1.0/clockrate)            
    !time=second()        
    
  end subroutine cputim



  ! the following routines are called at start-up to calculate some
  ! useful information for the model, like array adresses and 
  ! the lay-out of the defined regions.


  !
  ! determines children and ibeg,iend,jbeg,jend,lbeg,lend
  ! the subroutine is normally called once per run of TM5
  ! written by patrick berkvens and mike botchev, march-june 1999
  !

  subroutine determine_children_etc( status )

    use dims, only : nregions
    use dims, only : parent, children
    use dims, only : im, jm, lm
    use dims, only : xbeg, ybeg, zbeg
    use dims, only : xend, yend, zend
    use dims, only : xref, yref, zref
    use dims, only : ibeg, jbeg, lbeg
    use dims, only : iend, jend, lend
    use dims, only : isr, jsr
    use dims, only : ier, jer
    use dims, only : dx, dy, dz
    use dims, only : touch_np, touch_sp
    use dims, only : xcyc
    use dims, only : okdebug
    
    ! --- in/out ---------------------------------
    
    integer, intent(out)      ::  status
    
    ! --- const --------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/determine_children_etc'
 
    ! --- local ---------------------------------
    
    integer :: icells, jcells, lcells, jj, region
    integer :: xref_, yref_, zref_
    
    ! --- begin ---------------------------------

    if ( okdebug ) print *,'determine_children_etc: '

    children = 0
    do region=2,nregions
       ! increase children's counter of the parent:
       children(parent(region),0) = children(parent(region),0) + 1
       jj = children(parent(region),0)
       ! this child is region:
       children(parent(region),jj) = region
    end do
    ! for example:
    !children(0,0) = 1  ! means: there is one child in the region 0
    !children(0,1) = 1  ! means: this child is region 1
    !children(1,0) = 1  ! means: there is one child in region 1
    !children(1,1) = 2  ! means: this child is region 1
    !children(2,0) = 0  ! no children in region 2 

    do region=2,nregions
       ibeg(region) = nint((xbeg(region)-xbeg(parent(region)))    &
            /         (dx/xref(parent(region)))) + 1
       jbeg(region) = nint((ybeg(region)-ybeg(parent(region)))    &
            /         (dy/yref(parent(region)))) + 1
       lbeg(region) = nint((zbeg(region)-zbeg(parent(region)))    &
            /         (dz/zref(parent(region)))) + 1

       iend(region) = nint((xend(region)-xbeg(parent(region)))  &
            /         (dx/xref(parent(region)))) 
       jend(region) = nint((yend(region)-ybeg(parent(region)))  &
            /         (dy/yref(parent(region))))
       lend(region) = nint((zend(region)-zbeg(parent(region)))  &
            /         (dz/zref(parent(region))))

       ! test 
       xref_ = xref(region)/xref(parent(region))
       yref_ = yref(region)/yref(parent(region))
       zref_ = zref(region)/zref(parent(region))

       ! compute icells - how many cells of parent in x-direction a region occupies
       if (iend(region)>ibeg(region)) then
          icells = iend(region)-ibeg(region)+1
       else
          icells = im(parent(region)) - ibeg(region) + 1 + iend(region)
       end if
       if (icells*xref_/=im(region)) then
         write (gol,*) 'check xbeg/xend, im and ibeg/iend for region ',region; call goErr
         write (gol,*) '  xbeg,xend,im,ibeg,iend:', &
               xbeg(region),xend(region),im(region),ibeg(region),iend(region); call goErr
         TRACEBACK; call goErr; status=1; return
       end if
       jcells = jend(region)-jbeg(region)+1
       if (jcells*yref_/=jm(region)) then
          write (gol,*) 'check ybeg/yend, jm and jbeg/jend for region ',region; call goErr
          write (gol,*) '  ybeg,yend,jm,jbeg,jend:', &
               ybeg(region),yend(region),jm(region),jbeg(region),jend(region); call goErr
          TRACEBACK; call goErr; status=1; return
       end if
       lcells = lend(region)-lbeg(region)+1
       if (lcells*zref_/=lm(region)) then
          write (gol,*) 'check zbeg/zend, lm and lbeg/lend for region ',region; call goErr
          write (gol,*) '  zbeg,zend,lm,lbeg,lend:', &
               zbeg(region),zend(region),lm(region),lbeg(region),lend(region); call goErr
          TRACEBACK; call goErr; status=1; return
       end if
       
    end do
    if ( okdebug ) then
       write (gol,*) ' ibeg(2:nregions): ',ibeg; call goPr
       write (gol,*) ' iend(2:nregions): ',iend; call goPr
       write (gol,*) ' jbeg(2:nregions): ',jbeg; call goPr
       write (gol,*) ' jend(2:nregions): ',jend; call goPr
       write (gol,*) ' lbeg(2:nregions): ',lbeg; call goPr
       write (gol,*) ' lend(2:nregions): ',lend; call goPr
       write (gol,*) ' parents: ',parent; call goPr
       do region=1,nregions
          write (gol,*) ' children(',region, &
               ',0:nregions)=',children(region,:); call goPr
       end do
    end if
    ! determine isr,ier,jsr,jer: the scope of the cells to be 
    ! processed per region.
    ! processes like vertical mixing, chemistry, sources/sinks 
    ! are processed on this scope:
    ! do i=isr(region),ier(region)
    ! do j=jsr(region),jer(region)

    isr(1) = 1  
    ier(1) = im(1)
    jsr(1) = 1
    jer(1) = jm(1)
    do region = 2, nregions
       xref_ = xref(region)/xref(parent(region))
       yref_ = yref(region)/yref(parent(region))
       isr(region) = xref_ + 1
       ier(region) = im(region)-xref_
       jsr(region) = yref_+1
       jer(region) = jm(region)-yref_
       if(xcyc(region)/=0) then
          isr(region) = 1
          ier(region) = im(region)
       end if
       if(touch_sp(region)==1) jsr(region) = 1
       if(touch_np(region)==1) jer(region) = jm(region)
    end do
    
    ! ok
    status = 0

  end subroutine determine_children_etc


  subroutine define_splitorderzoom
    !
    ! splitorderzoom is splitorder specified for each region 
    ! the subroutine is normally called once per run of TM5
    ! ATTN: subroutine was written based on assumption that a region can not
    ! have number less than number of its parent: region > parent(region)
    ! written by patrick berkvens and mike botchev, march-june 1999
    !
    use dims
    implicit none
    ! local
    integer :: i, j, j0, region, tref_
    !character,dimension(3):: reverse

    if ( okdebug ) then
      write (gol,*) 'define_splitorderzoom: '; call goPr
    end if
    splitorderzoom = ' '
    splitorderzoom(1,1:nsplitsteps) = splitorder
    if ( zoom_mode == 1 ) then
       ! zoom_mode==1 means:
       !  x y z       | z y x
       !  x y z z y x | z y x x y z
       !  ...
       do region=2,nregions
          tref_ = tref(region)/tref(parent(region))
          if ( (tref_ > 4 ) .or. ( tref_ == 3 ) ) then
             write (gol,*) 'define_splitorderzoom: ERROR region = ',region; call goErr
             write (gol,*) 'define_splitorderzoom: ', 'wrong value for tref(region): ',tref(region); call goErr
             write (gol,*) 'define_splitorderzoom: ', 'should be 1,2 or 4 times bigger than tref of its parent'; call goErr
             write (gol,*) 'define_splitorderzoom: ', 'use another value for parameter zoom_mode'; call goErr
             stop
          end if
          j0 = 1  ! step counter for parent
          j = 1   ! step counter for region
          do while ( splitorderzoom(parent(region),j0) /= ' ' )
             ! use step triple of parent tref_ times:
             do i=1,tref_
                if (mod(i,2)==1) then
                   ! copy step triple of parent:
                   !WP! changed 2 into n_oper
                   splitorderzoom(region,j:j+(n_operators-1)) = &
                        splitorderzoom(parent(region),j0:j0+(n_operators-1))
                else
                   ! step triple of parent in reverse order:
                   !    WP! changed 2 into n_op 
                   splitorderzoom(region,j:j+(n_operators-1)) = & 
                        reverse(splitorderzoom(parent(region), &
                        j0:j0+(n_operators-1)))
                end if
                !WP! changed 3 into n_ope...
                j = j + (n_operators)
             end do
             ! next step triple of parent  !WP! changed 3 into n_oper..
             j0 = j0 + (n_operators)
          end do
       end do

    else if ( zoom_mode == 2 ) then

       ! zoom_mode==2 means:
       !  x y z       | z y x
       !  x y z x y z | z y x z y x
       !  ...
       do region=2,nregions
          tref_ = tref(region)/tref(parent(region))
          j0 = 1  ! step counter for parent
          j = 1   ! step counter for region
          do while (splitorderzoom(parent(region),j0)/=' ')
             ! replicate step triple of parent tref_ times:
             do i=1,tref_
                splitorderzoom(region,j:j+(n_operators-1)) = &
                     splitorderzoom(parent(region),j0:j0+(n_operators-1))      
                j = j + (n_operators)
             end do
             ! next step  of parent
             j0 = j0 + (n_operators)
          end do
       end do
    else
       write (gol,*) 'define_splitorderzoom: wrong value for zoom_mode ',zoom_mode; call goErr
       stop
    end if

  contains

    function reverse(str) result(inv_str)
      !WP! changed to n_operators....
      character,dimension(n_operators):: str,inv_str
      integer ::i

      do i=1,(n_operators)
         inv_str(i)=str(n_operators+1-i)
      end do

    end function reverse

  end subroutine define_splitorderzoom


  ! ***
  

  subroutine print_zoom(region_dat)

    use dims

    ! --- in/out ----------------------------------
    
    type(region_data), intent(in)   :: region_dat(nregions)

    ! --- local ---------------------------------
    
    integer :: region, j
    integer ::  n

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

    write (gol,'("print_zoom")'); call goPr
    write (gol,'("print_zoom: nregions : ",i6)') nregions; call goPr
    write (gol,'("print_zoom: xbeg     : ",20i6)') xbeg; call goPr
    write (gol,'("print_zoom: xend     : ",20i6)') xend; call goPr
    write (gol,'("print_zoom: ybeg     : ",20i6)') ybeg; call goPr
    write (gol,'("print_zoom: yend     : ",20i6)') yend; call goPr
    write (gol,'("print_zoom: zbeg     : ",20i6)') zbeg; call goPr
    write (gol,'("print_zoom: zend     : ",20i6)') zend; call goPr
    write (gol,'("print_zoom: xref     : ",20i6)') xref; call goPr
    write (gol,'("print_zoom: yref     : ",20i6)') yref; call goPr
    write (gol,'("print_zoom: zref     : ",20i6)') zref; call goPr
    write (gol,'("print_zoom: im       : ",20i6)') im; call goPr
    write (gol,'("print_zoom: jm       : ",20i6)') jm; call goPr
    write (gol,'("print_zoom: lm       : ",20i6)') lm; call goPr
    write (gol,'("print_zoom: isr      : ",20i6)') isr; call goPr
    write (gol,'("print_zoom: ier      : ",20i6)') ier; call goPr
    write (gol,'("print_zoom: jsr      : ",20i6)') jsr; call goPr
    write (gol,'("print_zoom: jer      : ",20i6)') jer; call goPr

    write (gol,'("print_zoom: ibeg(2:nregions) : ",20i6)') ibeg; call goPr
    write (gol,'("print_zoom: iend(2:nregions) : ",20i6)') iend; call goPr
    write (gol,'("print_zoom: jbeg(2:nregions) : ",20i6)') jbeg; call goPr
    write (gol,'("print_zoom: jend(2:nregions) : ",20i6)') jend; call goPr
    write (gol,'("print_zoom: lbeg(2:nregions) : ",20i6)') lbeg; call goPr
    write (gol,'("print_zoom: lend(2:nregions) : ",20i6)') lend; call goPr

    write (gol,'("print_zoom: parents : ",20i6)') parent; call goPr
    do region=1,nregions
       write (gol,'("print_zoom: children(",i2,",0:nregions)=",20i3)') region, children(region,:); call goPr
    end do

    do region=1,nregions
       n = count(splitorderzoom(region,:)/=' ')
       write (gol,'("print_zoom: splitorderzoom(",i2,",:)=",100(a1))') &
            region, splitorderzoom(region,1:n); call goPr
    end do

    do region = 1,nregions
       write (gol,'("print_zoom: Zoomed ... region ",i3)') region ; call goPr
       do j=jm(region),1,-1
          write (gol,'("print_zoom: ",i3,a3,200i1)') j,'---',region_dat(region)%zoomed(:,j); call goPr
       end do
       write (gol,'("print_zoom: Edge   ... region ",i3)') region ; call goPr
       do j=jm(region),1,-1
          write (gol,'("print_zoom: ",i3,a3,200i1)') j,'---',region_dat(region)%edge(:,j); call goPr
       end do
    end do

  end subroutine print_zoom

 
  ! ***
  

  subroutine blank_zoom_region(region,region_dat)
    !
    ! determine which cells in a region are calculated in an overlying
    ! zoom region. Only the core of a zoom region is flagged. The
    ! interface cells are not flagged, since they are treated coarse
    ! ref: Krol et al, IAMAS paper, 2001
    ! 
    ! Example    region1:
    ! Zoomed  I2 represent the core-zoom region
    ! 1111111111111111111
    ! 1111111111111111111
    ! 1111222111111111111
    ! 1111222111111111111
    ! 1111222111111111111
    ! 1111222111111111111
    ! 1111111111111111111
    ! 1111111111111111111
    ! Edge     !indicates the interface cells
    ! 0000000000000000000
    ! 0003222300000000000
    ! 0001000100000000000
    ! 0001000100000000000
    ! 0001000100000000000
    ! 0001000100000000000
    ! 0003222300000000000
    ! 0000000000000000000
    !     3 means x+y, 2 means y, and 1 means only x interface cell.
    !
    use dims
    implicit none

    ! input/output
    integer, intent(in)               :: region
    type(region_data),intent(inout)   :: region_dat

    ! local
    integer                           :: i, j, n, n1, xref_, yref_
    integer,dimension(:,:),pointer    :: zoomed, edge

    zoomed => region_dat%zoomed
    edge   => region_dat%edge

    zoomed = region
    edge = 0
    do n=1,children(region,0)  !loop over the children
       n1 = children(region,n)    !child region

       do i=1,im(region)
          if( (i <= ibeg(n1) .or. i >= iend(n1)) .and. (xcyc(n1) /= 1) ) cycle
          do j=1,jm(region)
             if( (j <= jbeg(n1)) .and. (touch_sp(n1)/=1) ) cycle
             if( (j >= jend(n1)) .and. (touch_np(n1)/=1) ) cycle
             zoomed(i,j) = n1
          end do
       end do

       if ( xcyc(n1) /= 1 ) then
          edge(ibeg(n1),jbeg(n1):jend(n1)) = 1
          edge(iend(n1),jbeg(n1):jend(n1)) = 1
       end if
       if ( touch_sp(n1) /= 1 ) then
          edge(ibeg(n1):iend(n1),jbeg(n1)) = &
               edge(ibeg(n1):iend(n1),jbeg(n1)) + 2
       end if
       if ( touch_np(n1) /= 1 ) then
          edge(ibeg(n1):iend(n1),jend(n1)) = &
               edge(ibeg(n1):iend(n1),jend(n1)) + 2
       end if

    end do

    ! put the edge of the region itself to -1.....
    if ( parent(region) /= 0 ) then
       xref_ = xref(region)/xref(parent(region))
       yref_ = yref(region)/yref(parent(region))
       if(xcyc(region)/=1) then
          edge(1:xref_,:) = -1
          edge(im(region)-xref_+1:im(region),:) = -1
       end if
       if(touch_sp(region)/=1) edge(:,1:yref_) = -1
       if(touch_np(region)/=1) edge(:,jm(region)-yref_+1:jm(region)) = -1
    end if

    nullify(zoomed)
    nullify(edge)

  end subroutine blank_zoom_region


end module initexit
