ecearth3/sources/tm5mp/proj/bacchus/initexit.F90

!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (rcode/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (rcode> 0) then; TRACEBACK; action; return; end if
#define IF_NOTOK_MPI(action) if (ierr/=MPI_SUCCESS) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:     INITEXIT
!
! !DESCRIPTION: contains routines to initialize/finalize the model
!\\
!\\
! !INTERFACE: 
!
MODULE INITEXIT
  !
  ! !USES:
  !
  use GO, only : gol, goPr, goErr, goLabel
  use dims  ! WARNING: has var 'status'
  
  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: exitus
  public :: start
  public :: control_init
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter ::  mname = 'initexit'
  !
  ! !REVISION HISTORY: 
  !    6 Nov 2012 - Ph. Le Sager - new read_control subroutine
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------
  
CONTAINS
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DEFAULT_CNTL
  !
  ! !DESCRIPTION: provide default values of control variables
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DEFAULT_CNTL ( rcode )
    !
    ! !USES:
    !
    use datetime,    only : tau2date
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) ::  rcode
    !
    ! !REVISION HISTORY: 
    !    mh, 27-jun-1989 - 26-sep-1992
    !    mk, 21-dec-2002
    !    6 Nov 2012 - Ph Le Sager -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/default_cntl'
    integer :: i,k,n

    ! set calendar type
    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 ) then
       rcode=1
       write(gol,*)' GMT start time should be multiple of 3'; call goErr
       TRACEBACK; return
    end if
    
    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
    okdebug=.false.
    revert = 1

  END SUBROUTINE DEFAULT_CNTL
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EXITUS
  !
  ! !DESCRIPTION: terminate a model run
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EXITUS( rcode )
    !
    ! !USES:
    !
    use global_data  , only : free_fields, outdir
    use io_save,       only : write_save_file
    use restart,       only : rs_write
    use datetime,      only : tstamp
    use advectm_cfl,   only : done_cfl
#ifdef with_budgets
    use budget_global, only : done_budget_global
#endif
#ifdef with_ecearth_optics
    use ecearth_optics, only : ECEarth_Optics_Done
#endif
    use Partools,      only : isRoot
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) ::  rcode
    !
    ! !REVISION HISTORY: 
    !    6 Nov 2012 - Ph Le Sager -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer   ::  region
    real      ::  cpu3

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

    ! Save the model state. This routine is quite memory-consuming, and should
    ! not be called if not used.
    if(.not.rs_write) then 
       call Write_save_file( 'successful completion of run',trim(outdir)//'/save', rcode )
       IF_NOTOK_RETURN(rcode=1)
    endif

#ifdef with_budgets
    ! save budgets, print summary
    call done_budget_global ( rcode )
    IF_NOTOK_RETURN(rcode=1)
#endif

#ifdef with_ecearth_optics
    call ECEarth_Optics_Done( rcode )
    IF_NOTOK_RETURN(rcode=1) 
#endif

    ! free memory
    call free_fields

    if ( isRoot ) then
       write (gol,'(" ")'); call goPr
       write (gol,'("CFL info from advection:")'); call goPr
       write (gol,'(a,i4,f10.4)') '  x: nloop_max, xi', nloop_max(1,1), xi(1,1); call goPr
       write (gol,'(a,i4,f10.4)') '  y: nloop_max, xi', nloop_max(1,2), xi(1,2); call goPr
       write (gol,'(a,i4,f10.4)') '  z: nloop_max, xi', nloop_max(1,3), xi(1,3); call goPr
    end if

    ! cfl finished
    call Done_CFL

    IF ( isRoot ) THEN
       write (gol,'(1x)'); call goPr
       write (gol,'("program has terminated normally.")'); call goPr
       call cputim(cpu3)
       cpu3 = cpu3-cpu0
       write (gol,'(a," > number of timesteps   :",i8)')    rname, nstep      ; call goPr
       write (gol,'(a," > time-loop runtime [s] :",f16.2)') rname, cpu3       ; call goPr
       write (gol,'(a," > runtime/timesteps [s] :",f16.8)') rname, cpu3/nstep ; call goPr       
       write (gol,'(1x)'); call goPr
    END IF

  END SUBROUTINE EXITUS
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CONTROL_INIT
  !
  ! !DESCRIPTION: set control variables, either to default or read from rc file.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CONTROL_INIT( rcode )
    !
    ! !USES:
    !
    use GO,            only : TrcFile, Init, Done, ReadRc
    use GO,            only : TDate, NewDate, AnyDate
    use GO,            only : operator(+), operator(-)
    use GO,            only : goTranslate
    use datetime,      only : date2tau, tau2date, julday
    use global_data,   only : rcfile, inputdir, outdir
    use global_data,   only : fcmode, tfcday0
    use partools
#ifdef oasis3
    use tm5_prism,     only : PRISM_start_date
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: rcode
    !
    ! !REVISION HISTORY: 
    !    6 Nov 2012 - Ph Le Sager - v0
    !
    ! !REMARKS:
    !  - this is code taken off 'start' in order to read control parameters
    !    before 'start' is called, so they are available for processes inits.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
    character(len=*), parameter ::  rname = mname//'/control_init'

    character(len=32)   ::  stime
    type(TrcFile)       ::  rcF
    integer             ::  ccyy, mm, dd
    
    ! -------------------- begin --------------------
    
    call default_cntl( rcode )
    IF_NOTOK_RETURN(rcode=1)

    onROOT : if ( isRoot ) then

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

       ! forecast series ?
       call ReadRc( rcF, 'time.fc', fcmode, rcode, default=.false. )
       IF_ERROR_RETURN(rcode=1)

       ! read forecast day 0 ?
       if ( fcmode ) then

          ! read day: yyyy-mm-dd
          call ReadRc( rcF, 'time.fc.day0' , stime, rcode )
          IF_NOTOK_RETURN(rcode=1)
          call goTranslate( stime, '/-', ' ', rcode )
          IF_NOTOK_RETURN(rcode=1)
          read (stime,*,iostat=rcode) ccyy, mm, dd
          if ( rcode /= 0 ) then
             write (gol,'("reading ccyy mm dd from : ",a)') trim(stime); call goErr
             TRACEBACK; call goErr; rcode=1; return
          end if
          ! store day:
          tfcday0 = NewDate(year=ccyy,month=mm,day=dd)

       else
          ! dummy:
          tfcday0 = AnyDate()
       end if

       call ReadRc(rcF, 'time.ndyn_max', ndyn, rcode )
       IF_NOTOK_RETURN(rcode=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, rcode )
       IF_NOTOK_RETURN(rcode=1)

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

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

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

       ! name of output directory:
       call ReadRc( rcF, 'output.dir', outdir, rcode )
       IF_NOTOK_RETURN(rcode=1)

       ! Start time of the entire simulation (i.e. the first day of the first leg)
       call ReadRc( rcF, 'timerange.start' , stime, rcode )
       IF_NOTOK_RETURN(rcode=1)
       call goTranslate( stime, '/-:', ' ', rcode )
       IF_NOTOK_RETURN(rcode=1)
       read (stime,*,iostat=rcode) sdate_simulation
       if ( rcode /= 0 ) then
          write (gol,'("could not read start time from : ",a)') trim(stime); call goErr
          TRACEBACK; rcode=1; return
       end if

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

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

       ! 'target' time?
       idatet = idatee

       ! close:
       call Done( rcF, rcode )
       IF_NOTOK_RETURN(rcode=1)

    end if onROOT

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

    ! Init time/calendar
    !
    ! 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 date2tau(idatee,itaue)
    call date2tau(idatet,itaut)

    ! set time flags
    newyr=.true.
    newhour(:)=.true.
    newday=.true.
    newmonth=.true.
    newsrun=.true.

    ! step counter
    nstep0=0
    nstep=0

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

    ! remains from zooming capabilities - set for region 1 if still needed
    children = 0
    isr(1) = 1  
    ier(1) = im(1)
    jsr(1) = 1
    jer(1) = jm(1)
    splitorderzoom = ' '
    splitorderzoom(1,1:nsplitsteps) = splitorder
    
    rcode=0
    
  END SUBROUTINE CONTROL_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    START
  !
  ! !DESCRIPTION: initialization of a model run or its continuation
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine start( tread1, tread2, rcode )
    !
    ! !USES:
    !
    use GO,            only : TrcFile, Init, Done, ReadRc
    use GO,            only : TDate, NewDate, IncrDate
    use GO,            only : operator(+), operator(-), rTotal
    use global_data,   only : rcfile
    use tracer_data,   only : tracer_print, init_short
    use Meteo,         only : Meteo_Setup_Mass, Meteo_Setup_Other
#ifndef __GFORTRAN__
    use tracer_data,   only : init_non_zero
#endif
#ifdef with_budgets
    use budget_global, only : Init_budget_global
#endif
    use advectm_cfl,   only : Init_CFL
#ifndef without_advection
    use advectm_cfl,   only : Check_CFL
#endif
#ifdef with_ecearth_optics
    use ecearth_optics, only : ECEarth_Optics_Init
#endif
    ! to fill tracers: 
    use user_input,   only : user_input_start
    use restart,      only : Restart_Read,      Restart_Write
    use io_save,      only : read_save_file_30, read_save_file
    use io_save,      only : readhdfmmr,        read_mmix
    !
    ! !OUTPUT PARAMETERS:
    !
    type(TDate), intent(out)    ::  tread1, tread2
    integer,     intent(out)    ::  rcode
    !
    ! !REVISION HISTORY: 
    !    6 Nov 2012 - Ph Le Sager - took off reading of control parameters
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer             ::  n, region
    character(len=256)  ::  fname, fdir
    character(len=32 )  ::  key
    type(TrcFile)       ::  rcF
    type(TDate)         ::  tdyn, tr(2)
    real                ::  t1, t2

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

    call goLabel()

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

    ! initialise data for CFL routines
    CALL init_cfl

    write (gol,'("  ",a,":   initial meteo (pressure, air mass), optionally from restart file ...")') 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, rcode, isfirst=.true. )
    IF_NOTOK_RETURN(rcode=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
    ! occurs
    call Check_CFL( tread1, tread2, n, rcode )
    IF_NOTOK_RETURN(rcode=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, rcode )
    !if (rcode/=0) call escape_tm('ERROR in tracer') 
    ! setup (interpolate?) other meteo:    
    call Meteo_Setup_Other( tr(1), tr(2), rcode, isfirst=.true. )
    IF_NOTOK_RETURN(rcode=1)


    !
    ! ** INIT TRACERS ****************************
    !     
    write (gol,'("  ",a,":   init tracer fields (istart=",i2,")...")') rname, istart; call goPr

    IF (revert == 1) THEN

       select case (istart)
       case(1)
          !
          ! initial tracer fields are set to zero
          ! nothing to do, since TRACER_INIT already set them to 0.

#ifndef __GFORTRAN__
       case ( 2 )
          !
          ! initial tracer fields with a very small non-zero values
          !
          call INIT_NON_ZERO
#endif

       case ( 30 )
          !
          ! Read "save" files with MIXING RATIO and (option) SLOPES, SECOND
          ! MOMENTS, of TRANSPORTED TRACERS.
          ! 
          ! File are defined with fully qualified name in the rc-file with
          ! the following format:
          !
          !   start.30.<region-name> : <fully-qualified-filename>
          !
          ! Example:
          ! 
          !   start.30.glb600x400    : /data/save_files/mystuff_glb6x4.hdf
          !
          call Init( rcF, rcfile, rcode )
          IF_NOTOK_RETURN(rcode=1)

          ! loop over regions:
          do region = 1, nregions

             write (key,'("start.30.",a)') trim(region_name(1))
             call ReadRc( rcF, key, fname, rcode, default='file_name_empty' )
             IF_NOTOK_RETURN (rcode = 1)
             write (gol,*) 'Using save file names from rc-file start.30.* values'; call goPr

             CALL Read_save_file_30( region, fname, rcode )
             IF_NOTOK_RETURN (rcode = 1)

          end do

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

       case ( 31 )
          !
          ! Read mass of both transpoted *AND* short lived species from so
          ! -called "save file". No slopes, but regridding to model resolution
          ! is available, since grid information is also read.
          !
          call Init( rcF, rcfile, rcode )
          IF_NOTOK_RETURN(rcode=1)


          do region=1,nregions

             ! name of save file
             call ReadRc( rcF, 'start.31.'//trim(region_name(region)), fname, rcode )
             IF_NOTOK_RETURN(rcode=1)

             ! read all tracers
             call Read_save_file( region, fname, rcode )
             IF_NOTOK_RETURN(rcode=1)

             ! overwrite with TM4 fields ?
             call ReadRc( rcF, 'start.31.'//trim(region_name(region))//'.TM4', fname, rcode, 'none' )
             IF_ERROR_RETURN(rcode=1)

             ! key found ? then read:
             if ( fname /= 'none' ) then
                call Read_save_file( region, fname, rcode, tm4=.true. )
                IF_NOTOK_RETURN(rcode=1)
             end if

          end do

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

       case ( 32, 33 )
          !
          ! 32 = read from restart file: tracers mass only
          ! 33 = read from restart file: tracers *AND* air masses
          !
          ! Note that for 33, surface pressure and air mass (both available in
          ! the restart file) are also read in Meteo_Setup_Mass above. But not
          ! with 32.
          !
          call cpu_time(t1)

          call Restart_Read( rcode, tracer_mass=.true., air_mass=.true.)
          IF_NOTOK_RETURN(rcode=1)

          call cpu_time(t2)
          write (gol,*) "   time to read restart [s]: ", t2-t1 ; call goPr

       case ( 4 )
          !
          ! Initial tracer fields are obtained from a "saveoldfile" : see
          ! io_save for details. 
          !
          call Init( rcF, rcfile, rcode )
          IF_NOTOK_RETURN(rcode=1)

          call ReadRc( rcF, 'start.4.'//trim(region_name(1)), fname, rcode )
          IF_NOTOK_RETURN(rcode=1)
          
          call readhdfmmr( 1, fname, rcode )
          IF_NOTOK_RETURN(rcode=1)

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

       case ( 5 )
          !
          ! Transported tracer fields are obtained from a mmix file.
          ! Slopes, if used, are set to 0.
          !
          ! 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 (tiny(1.))
          !
          call Init( rcF, rcfile, rcode )
          IF_NOTOK_RETURN(rcode=1)

          call ReadRc( rcF, 'start.5.'//trim(region_name(1)), fname, rcode )
          IF_NOTOK_RETURN(rcode=1)

          call READ_MMIX(1,fname, rcode)
          IF_NOTOK_RETURN(rcode=1)

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

       case ( 9 )
          !
          ! USER DEFINED
          !                  
          call user_input_start( rcode )
          IF_NOTOK_RETURN(rcode=1)

       case default

          write (gol,'("unsupported istart : ",i2)') istart; call goErr
          TRACEBACK; rcode=1; return

       end select
       
       ! Ensure that non-tranported tracers are initialized
       if((istart==4.or.istart==5) .and.newsrun) then
          do n=1,nregions
             call init_short(n)
          enddo
       endif
       
    end if   ! forward run

    if ( okdebug ) then
       do region=1,nregions
          call tracer_print(region, "read init", rcode)
          IF_NOTOK_RETURN(rcode=1)
       end do
    end if
        
#ifdef with_budgets
    write (gol,'("  ",a,":   init budgets ...")') rname; call goPr
    call init_budget_global ( rcode )
    IF_NOTOK_RETURN(rcode=1)
#endif

#ifdef with_ecearth_optics
    write (gol,'("  ",a,":   init ecearth optics ...")') rname; call goPr
    Call ECEarth_Optics_Init( rcode )
    IF_NOTOK_RETURN(rcode=1)
#endif

    call cputim(cpu0) ; call goLabel() ; rcode=0

  END SUBROUTINE START
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CPUTIM
  !
  ! !DESCRIPTION: returns current value of the processor clock in seconds.
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine cputim(time )
    !
    ! !OUTPUT PARAMETERS:
    !
    real,intent(out) :: time
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer :: clock, clockrate

    call system_clock(clock, clockrate)    
    time = clock*(1.0/clockrate)            

  end subroutine cputim
  !EOC
  
END MODULE INITEXIT